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Introduction 

This  document  is  a  final  report  presenting  the  results  of  a  four  year  grant.  This 
introduction  includes  a  brief  description  of  the  subject,  goals,  and  scope  of  the  research  project, 
along  with  a  description  of  the  relationship  of  this  work  to  that  of  previous  researchers. 
Following  this  description  is  a  milestone  chart  showing  the  focus  areas  and  a  timeline  of  the 
work  accomplished.  The  introductory  section  ends  with  a  brief  summary  of  the  final  program 
accomplishments  followed  by  a  more  detailed  review  of  the  results  of  the  first  three  years  of  the 
program.  The  main  body  of  the  report  is  a  detailed  description  of  the  accomplishments  during 
the  final  year  of  the  program.  This  is  followed  by  some  concluding  remarks  addressing  the 
significance  of  the  research  accomplishments  and  anticipated  future  uses  of  these  results. 

Program  Description  and  Goals 

Hyperthermia  has  been 
shown  to  be  an  efficacious  adjuvant 
therapy  in  the  treatment  of  recurrent 
breast  cancer.  [1,2,3]  One  of  the 
most  flexible  and  proven  methods  to 
induce  hyperthermia  is  ultrasound. 

In  order  to  optimally  utilize  such 
clinical  devices  it  is  essential  that 
the  operator  have  the  ability  to 
accurately  plan  a  treatment.  The 
clinical  hyperthermia  situation  is 
diagrammed  in  Figure  1.  A 
treatment  plan  consists  of  a  defined 
motion  path  and  applied  transducer 
power  that  will  raise  the  temperature 
of  the  tumor  to  a  desired  set  point 
for  a  specified  time  while 
minimizing  the  effect  on 
surrounding  tissue.  The  treatment 
plan  must  be  adjusted  for  tissue 
geometry,  absorption,  and  desired 
temperature  history  specific  to  each 
patient.  Planning  a  treatment  begins 
with  the  accurate  estimation  of  the 
absorbed  power  field  created  by  the 
ultrasound  transducer.  Thus,  the 

ability  to  accurately  estimate  this  field  is  critical  to  producing  a  proper  and  effective  treatment 
plan. 

It  is  the  goal  of  this  research  to  provide  improved  numerical  modeling  of  ultrasound 
propagation  through  breast  tissue  and  the  post-mastectomy  chest  wall  for  use  in  patient  treatment 
p  anning  of  hyperthermia  cancer  treatments.  With  this  improved  model  the  clinician  can  plan  an 
u  trasound  hyperthermia  treatment  so  as  to  produce  the  required  thermal  dose  while  minimizing 
patient  pain  due  to  excessive  temperature  or  ultrasound-tissue  interactions.  A  clinical  system  is 


Patient 


Tumor 


Membrane 


Ultrasound  Transducer 

Figure  1 :  Diagram  of  the  clinical  situation  for  a 
hyperthermia  treatment  of  the  breast  or  chest  wall 
area. 
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being  constructed  to  combine  model  results  with  experiments  in  an  inverse  technique  to  estimate 
the  absorbed  power  field  in  the  breast  and  chest  wall  during  ultrasound  hyperthermia.  This 
improved  planning  will  aid  the  clinicians  in  providing  better  hyperthermia  treatments,  which  will 
improve  treatment  response  rates. 

The  system  under  development  consists  of  the  following  subsystems: 

1 .  An  ultrasonic  B-mode  imager  based  clinical  data  acquisition  system  used  to  obtain 
patient  anatomy  and  measure  the  attenuation  of  ultrasound  by  the  body  tissue. 

2.  A  clinical  technique  using  an  ultrasound  treatment  system  to  directly  measure  absorbed 
power  at  specific  points. 

3.  An  improved  model  of  ultrasound  propagation  through  breast  tissue. 

4.  An  inverse  technique  that  integrates  experimental  measurements  and  modeling  results  to 
obtain  the  “best"  prediction  of  the  ultrasound  power  deposition  field. 

The  B-mode  diagnostic  ultrasound  imager  in  Subsystem  1  is  used  to  construct  a 
geometric  model  based  upon  each  patient's  anatomy  and  to  measure  the  ultrasound  attenuation 
coefficient  at  selected  areas  within  the  tissue  region.  In  Subsystem  2  the  clinical  treatment 
system  is  used  in  a  scanning  protocol  to  locate  thermocouples  imbedded  in  the  breast  tissue.  The 
response  of  these  thermocouples  to  the  scanned  ultrasound  gives  a  direct  measure  of  the  power 
deposition  at  that  point.  Subsystem  3  is  an  ultrasound  propagation  and  power  deposition  model. 
For  this  model  we  have  considered  2  options: 

a)  Hybrid  Green’s  Function  Model  -  Ultrasound  propagation  from  the  transducer  into  the 
breast  is  modeled  with  a  hybrid  of  Green's  function  solutions  used  in  the  homogenous 
water  region  and  a  finite  element  solution  of  the  acoustic  wave  equation  in  the 
inhomogeneous  breast  region. 

b)  Parabolic  Model  -  This  is  a  forward  marching  parabolic  model  that  is  quite  fast,  but  does 
not  take  into  account  reflected  energy  at  interfaces. 

Subsystem  4  is  an  iterative  inverse  technique  that  optimizes  the  parameters  of  the  power 
deposition  model  by  forcing  the  model  results  to  match  the  experimental  measurements.  The 
model  is  used  to  predict  the  power  deposition  at  the  thermocouple  locations,  the  error  between 
model  and  experiment  is  measured,  and  the  error  values  are  used  to  adjust  the  model  parameters 
for  the  next  iteration. 

This  improved  power  deposition  model  will  be  useful  in  the  development  of  more 
accurate  path  planning  for  hyperthermia  treatment.  As  such  it  will  provide  a  clinically  useful 
tool  that  can  serve  as  the  base  for  the  development  of  more  accurate  treatment  path  planning, 
more  accurate  control  techniques,  and  eventually,  higher  success  rates  in  hyperthermia  treatment 
with  reduced  patient  discomfort  during  treatment. 
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Task  Description  and  Milestones: 
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Outline  of  Accomplishments: 

1*‘  Year  Accomplishments: 

Subsystem  1 :  Geometry  Acquisition  and  Attenuation  Measurement 

•  Acquisition  of  ultrasonic  B-scan  instrument. 

•  Deyelopment  of  interface  hardware  to  isolate  analog  A-line  information. 

•  Acquisition  of  high  speed  A/D  system  and  deyelopment  of  control  software. 

•  Inyestigation  of  reflection  based  acoustic  attenuation  measurement  techniques. 

Subsystem  2:  Experimental  SAR  Measurement 

•  Initial  discussions  and  yerification  of  equipment. 

Subsystem  3:  Model  Deyelopment 

•  Deyelopment  and  testing  of  a  2D  finite  element  based  model. 

•  Identification  of  altematiye  solutions. 

Subsystem  4:  Inyerse  Technique 

•  Initial  discussions. 

2"**  Year  Accomplishments: 

Subsystem  1 :  Geometry  Acquisition  and  Attenuation  Measurement 

•  Deyelopment  of  Log  Spectral  Difference  software  for  measuring  the  acoustic 
attenuation  from  a  single  A-line. 

•  Deyelopment  and  testing  of  a  Force/Balance  test  apparatus  for  direct  measurement  of 
acoustic  attenuation. 

Subsystem  2:  Experimental  SAR  Measurement 

•  No  work  accomplished. 

Subsystem  3:  Model  Deyelopment 

•  Extension  of  the  2D  FE  model  to  three  dimensions.  Combination  of  the  FE  model 
with  an  integral  method  model  to  improve  the  speed  of  the  model. 

•  Development  of  a  3D  parabolic  model  as  a  much  faster,  possibly  less  accurate 
alternative. 

•  Development  of  a  3D  Green’s  Function  model  for  use  as  the  “gold  standard”  with 
which  the  other  models  will  be  compared. 

•  Initial  model  comparisons. 

Subsystem  4:  Inverse  Technique 

•  No  progress. 

Year  Accomplishments: 

Subsystem  1 :  Geometry  Acquisition  and  Attenuation  Measmement 

•  Laboratory  hardware  was  developed  to  provide  3D  scanning  capability  for  the 
medical  B-scan  imager.  Clinical  hardware  is  now  essentially  complete. 

•  Software  to  produce  a  3D,  tissue  attenuation  map  from  the  A-mode  signals  was 
completed  and  is  undergoing  verification  testing. 
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•  Software  and  hardware  for  clinical  patient  geometry  acquisition  is  approximately 
50%  complete. 

Subsystem  2:  Experimental  SAR  Measurement 

•  Thermocouple  based  beam  plots  have  been  performed  to  calibrate  the  ultrasound 
power  field  in  water  for  specific  transducers. 

•  Experiments  have  been  conducted  in  water  and  are  beginning  using  cast  agar 
phantoms. 

•  Preliminary  in-vivo  experiments  were  conducted  on  a  dog. 

•  A  finite  difference  model  has  been  constructed  to  help  interpret  the  experimental 
measurements. 

Subsystems:  Model  Development 

•  The  parabolic  model  has  been  selected  for  use  in  regions  of  low  acoustic  contrast.  In 
these  regions  the  model  is  quite  fast  and  loses  little  accuracy  when  compared  to  the 
Green’s  Function  Model.  The  Green’s  Function  Model  can  be  used  in  high  contrast 
areas  (close  to  ribs,  etc.) 

•  User  interface  developed  allowing  specification  of  the  volume  modeled,  the  acoustic 
power  input  as  a  2D  complex  pressure  function,  and  model  parameters;  density, 
sound  velocity,  and  absorption  at  each  node. 

Subsystem  4:  Inverse  Technique 

•  Basic  algorithm  development.  Optimization  runs  on  simple  geometries. 

4'*'  Year  Accomplishments: 

Subsystem  1 :  Geometry  Acquisition  and  Attenuation  Measurement 

•  Verification  tests  for  clinical  system  for  3D  attenuation  mapping. 

•  Software  and  hardware  for  clinical  patient  geometry  acquisition  completed. 

•  Verification  tests  completed. 

Subsystem  2:  Experimental  SAR  Measurement 

•  Probe  orientation  effects  on  the  SAR  measurements  were  investigated 
experimentally. 

•  A  series  of  new  probes  were  designed,  constructed  and  tested. 

•  A  final  design  was  proposed  for  future  work. 

Subsystems:  Model  Development 

•  The  model  was  completed  in  the  3'"*  year.  No  work  was  accomplished  on  this 
subsystem  during  the  4*  year. 

Subsystem  4:  Inverse  Technique 

•  Subsystem  integration  completed. 

•  Technique  used  successfully  to  optimize  the  input  power  model. 

•  Given  the  level  of  scatter  in  the  experimental  SAR  measurements,  the  optimization 
model  was  not  used  to  adjust  nodal  attenuation  values. 
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Program  Summary 

The  four  year  research  effort  has  resulted  in  an  improved  power  deposition  model  for 
hyperthermia  treatment  planning.  The  model  is  significantly  faster  than  previous  power 
deposition  models  used  for  treatment  planmng.  It  is  more  flexible,  accepting  tissue  geometry 
and  attenuation  values  measured  clinically.  It  uses  an  inverse  optimization  scheme  to  improve 
die  fit  between  the  model  parameters  and  clinical  measurements.  This  optimization  has  been 
implemented  for  timing  the  power  input  model  and  can  potentially  be  implemented  for  tuning 
tissue  property  data.  Extensive  testing  has  shown  that  existing  thermocouple  probes  are  not  well 
suited  to  measure  SAR.  An  improved  thermocouple  has  been  designed,  but  construction,  testing, 
and  implementation  will  be  left  to  future  researchers.  This  system  is  now  being  integrated  into 
the  hyperthermia  treatment  facility  at  the  University  of  Utah  Hospital. 

During  the  first  two  years  of  the  research  program,  a  parabolic  numerical  model  was 
developed  that  is  dramatically  faster  than  finite  element  models  or  models  based  on  integral 
equations.  The  model  is  capable  of  processing  a  model  with  10-^6  nodes  in  approximately  2 
minutes  (using  a  Sun  Sparc  20  workstation).  The  model  handles  free  space  boundary  conditions 
much  more  conveniently  than  the  finite  element  model.  The  model  accurately  handles  the 
forward  propagation  of  ultrasonic  energy,  but  does  not  account  for  reflected  energy.  This  is  not 
a  significant  limitation  since  the  low  acoustic  contrast  occurring  in  tissue  does  not  give  rise  to 
large  reflections.  Comparisons  to  an  integral  model  showed  a  good  match  except  when  bone 
was  present.  This  is  considered  acceptable  since  treatment  requiring  significant  insonation  near 
bones  is  unlikely  due  to  patient  pain.  These  results  are  being  submitted  for  journal  publication. 

The  system  includes  an  ultrasonic  B-scan  based  patient  geometry  acquisition  and 
attenuation  measurement  system  that  allows  the  model  nodal  parameters  to  be  conveniently 
adjusted  to  match  the  tissue  geometry  and  attenuation  of  each  patient.  This  system  requires  a 
pre-treatment  scanning  session  where  an  operator  acquires  the  patient  geometry,  a  technician 
identifies  the  significant  tissue  regions,  and  the  attenuation  calculations  are  performed.  The 

system  has  been  shown  to  identify  tissue  geometry  to  an  accuracy  of +/- 1  mm  and  attenuation 
values  to  within  10%. 

An  inverse  optimization  technique  was  developed  to  adjust  model  parameters  to 
minimize  the  error  between  model  predicted  SAR  values  and  experimental  measurements.  The 
technique  was  developed  successfully,  however,  accurate  experimental  SAR  measurements  have 
not  been  achieved.  The  model  has  been  used  successfully  to  generate  an  accurate  model  for 
input  ultrasound  energy.  Here  the  inverse  technique  was  used  to  adjust  the  surface  velocity 
profile  of  a  model  of  the  ultrasound  transducer  to  provide  an  optimal  match  to  measurements 
made  with  a  hydrophone  probe.  Experimental  results  showed  that  when  a  transducer  model 
using  a  piston  like  velocity  profile  was  replaced  with  the  same  model  using  an  optimized  profile, 
the  error  between  experimental  hydrophone  measurements  and  the  model  prediction  was 
reduced  by  over  50%.  These  results  are  being  submitted  for  journal  publication. 

This  research  also  included  development  of  a  method  for  experimentally  measuring  the 
SAR  at  various  points  in  the  treatment  region  from  the  rate  of  temperature  rise  in  the  plastic 
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sheath  of  the  clinical  thermocouple  probes  currently  used  for  hyperthermia  treatment.  Results  of 
the  experiments  show  that  the  probe  measurement  is  dependent  on  both  the  angle  of  incidence  of 
the  ultrasound  wave  and  the  angle  of  rotation  of  the  probe.  The  angle  of  incidence  variation  has 
been  documented  in  the  literature,  but  the  angle  of  rotation  variation  was  a  surprise.  This 
variation  can  be  attributed  to  angularly  dependent  contact  between  the  thermocouple  and  the 
probe  sheath.  A  new  probe  design  has  been  designed,  but  has  not  yet  been  constructed  or  tested. 
Examining  the  experimental  results,  it  is  anticipated  that  this  technique  will  produce 
measurement  errors  of  not  less  than  +/-  5%. 

The  results  of  this  project  represent  a  significant  step  forward  in  accuracy  for  a  practical, 
clinical  ultrasound  power  deposition  model  for  hyperthermia  treatment  planning.  The  system  is 
currently  being  implemented  on  the  hyperthermia  treatment  facility  at  the  University  of  Utah 
Hospital  and  clinical  trials  are  anticipated  in  future  work. 


First  Year  Summary 

The  first  year  efforts  centered  on  initial  development  of  the  ultrasound  propagation 
models  and  acquisition  and  initial  development  of  the  diagnostic  ultrasound  equipment.  A  brief 
discussion  of  these  efforts  is  included  below.  For  more  detailed  information  refer  to  the  first 
year  progress  report. 

Initial  modeling  efforts  focused  on  two  potential  methods:  a  Green’s  Function  Method 
where  the  transducer  face  is  modeled  as  a  distribution  of  point  sources,  and  the  Finite  Element 
Method  where  FE  techniques  are  used  to  solve  the  acoustic  wave  equation.  A  two  dimensional 
FE  code  was  written,  and  verification  runs  were  under  way  at  the  end  of  the  first  year.  The 
results  of  these  runs  and  further  investigation  of  the  modeling  possibilities  led  to  the 
development  of  a  hybrid  model  that  combines  the  features  of  the  two  techniques  to  provide 
improved  speed  and  accuracy.  Details  of  this  model  are  discussed  in  the  second  year 
accomplishments  below. 

During  the  first  year  a  diagnostic  ultrasound  B-mode  imager  was  acquired.  This  imager 
was  modified  to  allow  access  to  the  raw  RF  waveform  (A-line)  for  each  scan  line  of  the  B-mode 
image.  A  schematic  of  the  resulting  system  is  shown  in  Figure  2.  The  ultrasound  unit  (Scanner 
250,  Pie  Medical  USA  3535  Route  66,  Neptune  NJ  07753)  was  modified  to  output  the  RF 
signals  along  with  sync  signals.  A  custom  hardware  selection  device  was  constructed  to  read 
these  signals  and  output  a  single  A-line.  A  PC  based  data  acquisition  system  was  acquired  that 
allows  digitization  of  the  A-line  signals.  Once  the  data  is  in  the  PC,  it  is  available  for  use  in 
software  techniques  to  estimate  the  ultrasonic  attenuation  field  in  the  tissue.  Along  with  the 
hardware  development,  software  was  developed  to  operate  the  hardware  and  perform  the 
attenuation  estimate.  [4] 

Second  Year  Summary 

Subsystem  1:  Ultrasound  Based  Anatomy  and  Attenuation  Measurement 
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Ultrasonic  Attenuation  Measurement 

The  power  deposited  at  a  specific  point  in  tissue  by  the  ultrasound  transducer  depends  on 
the  intensity  of  the  sound  field  at  that  point  and  the  absorption  coefficient  of  the  tissue  at  that 
point.  Clearly,  the  intensity  at  a  point  depends  on  the  initial  intensity  and  the  attenuation  of  the 
signal  in  the  tissue.  This  attenuation  is  due  to  scattering,  absorption,  and  geometric  attenuation. 
Previous  research  has  shown  that  the  attenuation  is  highly  variable  within  human  body  tissue, 
with  measurements  varying  significantly  for  different  tissue  types,  for  similar  tissue  on  different 
subjects,  and  even  for  the  same  tissue  on  the  same  patient  on  different  days.  [6-9]  Thus,  it  is 
important  to  measure  the  attenuation  of  the  ultrasonic  signal  in  vivo  on  the  patient  at  the 
treatment  site. 

Using  the  data  acquisition  system  shown  in  Fig.  2,  it  is  possible  to  estimate  the 
attenuation  of  the  ultrasound  at  different  points  in  the  patient  tissue  from  the  RF  signal  collected. 
This  is  accomplished  using  a  Log  Spectral  Difference  technique.  [7,8, 10]  Details  of  this 
technique  are  discussed  in  the  second  year  progress  report.  It  is  important  to  note  that  Log 
Spectral  Difference  is  one  of  several  possible  techniques  and  provides  an  estimate  of  the 
attenuation.  [5]  For  this  reason  it  is  important  to  test  the  accuracy  of  the  technique  with  other 
measurements. 

Tests  were  conducted  where  the  acoustic  attenuation  of  a  series  of  samples  were 
measured  using  both  log  spectral  difference  and  the  Force  Balance  technique.  Using  the  Force 
Balance  technique  the  force  generated  when  an  acoustic  wave  impinges  on  an  absorbing  material 
is  measured  and  related  to  the  power  of  the  wave.  Using  this  technique  it  is  possible  to  make  a 
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direct  measurement  of  the  attenuation  or  absorbed  power.  For  this  reason,  the  force  balance 
technique  is  used  as  a  standard  to  which  the  log  spectral  difference  measurements  are  compared. 
While  it  is  generally  considered  accurate  and  reliable,  the  force  balance  technique  is  not  suitable 
for  use  in  the  clinical  experiments  due  to  geometric  constraints. 

Using  the  techmques  described  above,  initial  tests  have  been  run  on  chicken  breast  to 
compare  the  results  of  the  two  techniques  discussed  above.  Using  the  Log  Spectral  Difference 
Technique,  the  normalized  attenuation  coefficient  for  the  chicken  was  measured  as  0. 1608 
Nepers  x  cm^  x  Using  the  Force  Balance  system  measured  values  ranged  from  0. 108 

to  0. 138  Nepers  x  crri^  x  MHz’^ .  These  values  are  of  the  same  magnitude  as  values  reported  in 
the  literature  (0.06  -  0. 13  Nepers  x  cm'^  x  1] 

Subsystem  2:  Absorbed  Power  Measurement 
No  work  accomplished  in  year  2. 

Subsystem  3:  Ultrasound  Propagation  Modeling 

During  the  second  year  two  competing  models  were  developed  and  compared.  The  first 
model  was  a  hybrid  numerical  model  using  integral  methods  to  solve  for  the  acoustic 
displacements  in  the  region  where  the  ultrasound  propagates  from  the  transducer  to  the  patient’s 
skin  and  using  a  finite  element  solution  of  the  wave  equation  for  propagation  in  the  body.  The 
model  predicts  the  displacements  and  stresses  propagating  through  the  tissue,  automatically 
accounting  for  reflections  and  transmission  at  all  interfaces.  This  model  was  expected  to  be 
quite  accurate.  A  second  model,  based  on  the  parabolic  equation  was  also  developed.  This 
model  accounts  for  forward  scattering  and  refraction,  but  ignores  reflections.  It  was  expected  to 
be  much  faster,  but  it  was  unclear  if  this  model  would  provide  sufficient  accuracy. 

Results:  Hybrid  Model 

The  Hybrid  Finite  Element  model  presented  two  difficulties; 

1 .  To  keep  the  computational  time  reasonable,  it  was  necessary  to  model  a  relatively 
small  volume.  It  was  not  possible  to  produce  boundary  conditions  at  the  edge  of  the 
volume  that  approximated  free  space  and  the  model  results  included  significant 
standing  waves  generated  at  the  model  boundaries.  Expanding  the  model  to  include 
the  entire  chest  cavity  would  solve  this  problem,  but  the  model  would  become 
unmanageably  large. 

2.  Initial  tests  were  performed  at  reduced  frequency  with  success.  When  frequencies 
corresponding  to  the  treatment  transducers  were  used,  the  finer  mesh  required 
resulted  in  extremely  high  computation  times. 
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Results:  Parabolic  Model 

The  parabolic  equation  method  is  an  efficient 
method  for  modeling  acoustic  wave  propagation  through 
low  contrast  acoustic  materials.  This  solution  technique  is 
somewhat  more  of  an  approximation  than  the  hybrid 
model,  but  because  of  the  geometry  used  in  the 
hyperthermia  experiment,  and  the  relatively  low  contrast  of 
breast  tissue,  it  is  possible  to  utilize  this  efficient 
approximation.  The  parabolic  model  was  found  to  be 
computationally  8  time  more  efficient  than  the  Green’s 
function  model  and  64  times  more  efficient  than  the  hybrid 
model. 


Figures  3  and  4  demonstrate  the  kind  of 
information  available  from  the  models.  In  Figure  3  we  see 
a  continuous  plane  wave  propagating  from  the  left, 
entering  the  body  through  a  vertical  skin  layer.  The  wave 
travels  through  an  elliptical  fatty  region  and  into  the 
spherical  tumor.  The  top  figure  shows  the  wave  energy 
density  that  tends  to  get  weaker  as  the  wave  moves  to  the 
right.  Note  that  the  tumor  absorbs  most  of  the  wave  energy 
resulting  in  a  “dark”  shadow  behind  the  tumor.  The 
bottom  figure  shows  the  energy  absorbed.  Note  that  in  the 
water  and  normal  tissue  the  absorbed  energy  is  low;  even 
when  the  wave  energy  is  high.  There  is  an  increase  in 
absorbed  energy  as  the  wave  passes  into  the  tumor.  This 
results  from  the  higher  absorption  coefficient  in  the  tumor. 
This  effect  is  shown  more  clearly  in  Fig.  4  which  is  a  line 
plot  of  the  absorbed  energy  along  a  central  axis  in  Fig.  3. 

Model  Accuracy 


Figure  3:  The  upper  image  is  the  incident 
energy  density.  The  lower  image 
shows  energy  absorbed  by  the  tissue. 


To  test  the  accuracy  of  the 
parabolic  method,  results  were  compared 
to  an  integral  equation  solution  of  the 
same  problem  that  is  representative  of  the 
best  models  currently  in  use.  The 
following  figures  demonstrate  the 
agreement  between  the  integral  equation 
method  and  the  parabolic  method.  Figure 
5  shows  the  geometry  used  for  this 
analysis.  Figure  6a  shows  a  color  mapped 
representation  of  a  slice  of  the  field.  Figure 
6b  is  an  axial  plot  of  the  same  data. 

Figure  7a  and  7b  show  the  results  using  a 
speed  of  sound  in  the  sphere  significantly 
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Figure  4:  Thermal  energy  deposition  along  a  central  axis 

predicted  by  the  parabolic  model. 
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less  than  and  significantly  greater  than  that  of 
water.  Given  the  differences  in  the  two  models, 
one  would  expect  error  to  increase  as  the 
mismatch  in  sound  velocity  between  the  tumor 
and  water  increases.  Thus,  these  two  conditions 
should  represent  worst  case  errors. 

The  excellent  match  between  the  two 
models  in  the  forward  scattered  diffraction 
patterns  shown  in  Fig.  7  and  the  poor  match  in  the 
reflected  dif&action  pattern  shown  in  Fig.  6b 
highlight  the  differences  between  the  two  models. 

The  excellent  overall  match  between  the  two 
models  in  these  worst  case  conditions,  coupled  with  the  speed  provided  by  the  parabolic  model 
make  it  an  excellent  candidate  for  our  application. 

Comparison  of  axial  field  I 
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Figure  6:  a)  Color  mapped  image  of  the  power  deposition  in  the  model  from  Fig.  5. 

b)  The  integral  equation  is  shown  in  solid,  and  shows  the  diffraction  pattern  due  to  reflections  from  the  hard 
sphere.  This  plot  shows  quantitatively  that  the  parabolic  method  does  not  take  into  account  the  reflections 
of  sound  waves.  Conditions  for  this  graph  correspond  to  those  in  Fig.  7b. 
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Figure  7:  a)  The  speed  of  sound  internal  to  the  sphere  was  1267.7  m/sec,  substantially  less  than  that  which  would  be 
encountered  in  the  breast,  thereby  making  the  impedance  contrast  quite  a  bit  larger  than  is  likely  to  be 
encountered  in  practice,  b)  Here,  the  speed  of  sound  internal  to  the  sphere  was  1677  m/sec,  which  is 
substantially  greater  than  any  tumor  that  would  be  found  in  the  breast.  In  each  case  the  solid  line  represents  the 
integral  model  and  the  dashed  line  represents  the  parabolic  model. 
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Subsystem  4:  Inverse  Optimization  Technique 
No  work  accomplished  during  the  2nd  year. 

Third  Year  Summary 

Subsystem  1;  Ultrasound  Based  Anatomy  and  Attenuation  Measurement 

•  Laboratory  hardware  was  developed  to  provide  3D  scanning  capability  for  the 
medical  B-scan  imager. 

•  Software  to  produce  a  3D,  tissue  attenuation  map  from  the  A-mode  signals  was 
completed. 

•  Software  and  hardware  for  clinical  patient  geometry  approximately  50%  complete. 


Patient  Geometry  Acquisition  and  Registration 

The  goal  here  is  to  be  able  to  place  a  patient  in  a  known,  repeatable  location  on  the 
treatment  station  and,  using  the  B-scan  imager,  identify  the  geometric  structure  of  the  tissue  and 
tumor  in  the  region  imder  treatment. 

Clinical  System  Description 

The  patient  positioning  system  consists  of  several  integrated  steps.  The  first  step  is  a 
repeatable  patient  positioning  system.  While  concepts  were  developed  during  this  year  no  final 
design  or  construction  were  completed. 

In  the  second 
phase  of  positioning  the 
patient,  a  series  of  B-scan 
ultrasound  images  of  the 
treatment  area  are 
collected.  These  images 
are  used  to  identify  the 
locations  of 

thermocouples,  interfaces 
between  different  tissue 
types,  and  other 
anatomical  features 
within  the  region  that 
will  be  affected  by  the 
hyperthermia  treatment. 

The  B-scan  images  are 
digitized  using  a  video 
capture  card.  A  custom 
software,  developed  as 
part  of  this  research. 


Figure  8:  The  software  under  development  acquires  a  series  of  B-scan  images, 
tiles  them  together  creating  a  larger  2D  image,  and  coordinates  the 
display  of  the  2D  images  to  represent  3D  geometry. 
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coordinates  multiple  B-scan  images  into  a  series  of  slices  representing  the  3D  patient  geometry. 
A  screen  capture  from  this  system  is  shown  in  Figure  8. 

The  doctor  or  technician  traces  features  of  the  B-scan  images  that  are  of  interest.  Using 
an  attenuation  map  developed  using  a  process  described  below  as  a  reference,  the  different 
regions  identified  by  the  tracings  are  assigned  attenuation  values.  The  software  processes  this 
information  and  produces  a  3D  grid  of  nodes  where  each  node  is  assigned  appropriate  values  to 
represent  the  patient  geometry  for  processing  by  the  math  model.  Figure  9  demonstrates  this 
process. 

The  left  most  image  contains  3  B-scans  that  have  been  registered  and  combined, 
producing  a  single  image.  In  the  center  image,  the  doctor  or  technician  has  mouse  sketched 
different  regions.  In  the  right  most  image,  these  regions  have  been  assigned  different  acoustic 
properties  based  on  apriori  knowledge  and  information  from  the  attenuation  mapping  system 
discussed  below. 

The  following  outline  describes  the  anticipated  clinical  procedure  using  the  developed 
system.  Depending  on  the  time  required  to  complete  the  inverse  optimi2ation  procedure,  the 

patient  may  either  lay  in  position  and  wait,  or  be  released  and  return  for  treatment  hours  or  days 
later. 


Figure  9:  Each  image  is  processed  by  a  doctor  or  technician.  The  process  involves  mouse  sketching  the  borders  of 
Identifiable  regions  and  assigning  each  region  uniform  acoustic  values  estimated  from  an  attenuation  map. 
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Three-dimensional  Ultrasonic  Attenuation  Field  Measurement 

As  mentioned  above,  the  goal  of  this  research  is  to  produce  a  model  that  will  accurately 
predict  the  power  deposited  at  all  points  in  a  region  of  tissue,  during  a  hyperthermia  treatment. 
This  will  allow  the  physician  to  develop  a  better  treatment  plan  that  will  maximize  the  positive 
effects  and  minimize  pain  and  other  adverse  effects.  In  order  for  the  inverse  model  to 
accomplish  this  efficiently,  the  model  must  begin  with  model  parameters  that  are  “in  the  ball 
park”  of  the  final  optimal  values.  The  least  well  known  of  these  parameters  is  the  acoustic 
attenuation  of  the  tissue.  Thus,  a  technique  was  developed  that  predicts  the  acoustic  attenuation 
field  within  the  tissue  based  on  data  fi-om  a  B-scan  imager.  The  third  year  efforts  extended  the 
previous  accomplishments  to  include  the  ability  to  produce  a  3D  attenuation  map  from  a  series  of 
B-scan  images. 

Background 

Measurement  of  ultrasonic  attenuation  from  reflected  signals  has  shown  promise  of  being 
a  valuable  diagnostic  tool[  12-20].  It  allows  physicians  to  diagnose  illnesses  that  cannot  be 
diagnosed  on  a  B-mode  ultrasound  image.  Because  of  this,  several  methods  of  calculating  the 
ultrasonic  attenuation  from  reflected  signals  have  been  developed  [14,17,19-22].  These  methods 
were  then  used  to  calculate  the  global  attenuation  of  a  structure.  The  next  step  in  attenuation 
calculation  was  to  produce  an  attenuation  image  that  could  be  used  in  conjunction  with  the  B- 
mode  image.  Walach  and  others  have  developed  techniques  that  relate  a  two-dimensional  B- 
mode  image  to  a  two-dimensional  attenuation  image  or  map  [12,13,23].  Seeing  both  images 
side-by-side,  the  physician  has  an  increased  diagnostic  capability. 

However,  the  usefulness  of  an  attenuation  map  is  not  limited  to  diagnostics.  A  3D 
attenuation  mapping  technique  has  been  developed  here  that  not  only  provides  the  initial 
estimates  for  the  power  deposition  model.  But  has  potential  to  ]»oduce  enhanced  resolution  2D 
attenuation  images. 

Method 

In  order  to  build  the  three-dimensional  attenuation  map,  a  B-mode  imaging  transducer 
was  placed  on  a  two-dunensional  linear  motion  machine.  The  transducer  is  part  of  an  ultrasonic 
imaging  system  (PIE  Scanner  250)  which  provides  access  to  the  reflected  ultrasound  signal  (A- 
line),  before  processing,  as  well  as  synchronization  signals  to  give  the  location  of  the  waveform 
along  the  width  of  the  transducer.  With  this  in  place,  software  was  developed  that  graniK-^t  the 
transducer  to  cover  the  entire  geometry  of  the  structure  being  modeled,  while  collecting  and 
organizing  the  data.  This  is  the  same  data  shown  in  Figures  8  and  9  collected  by  the  patient 
registration  system,  except  that  the  raw  A-lines  are  recorded. 

The  data  was  collected  with  an  eight-bit  analog  to  digital  converter,  and  was  then  written 
to  disk.  A  technique  has  been  developed  that  uses  several  different  gain  settings  for  each 
collection  to  mimic  a  higher-precision  conversion.  Using  this  technique,  each  A-line  was 
analyzed  at  the  highest  possible  gain  setting  where  saturation  did  not  occur.  The  focus  of  the 
transducer  was  set  so  that  the  structure  being  analyzed  was  in  the  far  field  of  the  transducer.  This 
insures  that  diffraction  effects  in  the  near  field  do  not  enter  into  the  calculations.  Once  the  data 
was  collected  it  was  then  analyzed  to  produce  the  attenuation  map. 


The  log  spectral  difference 
method  was  used  to  estimate  the  slope 
of  the  attenuation  to  frequency 
relationship.  (Lyons  and  Parker 
found  that  for  most  soft  tissues,  the 
attenuation  was  related  to  fi%quency 
to  the  power  of  n,  where  n  is  equal  to 
1 .0  - 1 .3  [25].  The  value  for  n  was  set 
to  1  in  this  research.)  This  method 
utilizes  the  ratio  of  the  power 
spectrum  from  a  shallow  and  deep 
segment  of  the  region  of  interest  (see 
Figure  10)  to  compute  the  attenuation 
slope  for  that  region.  The  specific 
steps  of  this  technique  have  been 
presented  in  the  2"^  year  progress 
report  as  well  as  previously  in  the 
literature  and  will  not  be  discussed 
here  [ 1 4, 1 5,22,26].  Our  application 
of  this  method  for  the  production  of  a  3D,  attenuation  map  will  be  discussed  next. 

As  mentioned  previously,  the  log  spectral  difference  technique  is  a  way  of  estimating  the 
slope  of  the  attenuation  for  a  specified  Region  Of  Interest  (ROI).  To  produce  die  attenuation 
map,  a  ROI  of  a  determined  size  was  moved  throughout  die  three-dimensional  data  space  until  an 
attenuation  map  was  produced.  The  region  of  interest  selected  was  1 .0  cm  wide  and  1 .7  cm  in 
depth  (see  Figure  10).  The  width  was  chosen  because  it  allowed  the  computation  to  be  averaged 
over  ten  reflected  signals.  The  depth  for  the  region  of  interest  came  fixim  a  formula  presented  by 
Kuc  [26]  and  the  specifications  of  the  transducer.  Once  the  ROI  size  was  selected,  it  was  moved 
throughout  the  data  space 
incrementally  in  1.0  mm  steps  in 
all  three  directions.  The 
attenuation-slope  value  for  each 
calculation  was  assigned  to  a 
point  in  the  center  of  the  ROI, 
thus  producing  an  attenuation 
map.  This  2D  ROI  is  similar  to 
that  used  by  previous 
researchers.  The  additional 
information  available  from  the 
3D  data  collection  used  here 
allows  a  3D  ROI  to  be 
investigated,  where  data  from  a 
3D  region  are  averaged  together 
to  improve  the  accuracy  and/or 
resolution  of  the  measurement. 

Software  has  been  developed  to 
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completely  automate  the  task  of  data  acquisition  and  analysis,  and  calibration  tests  are  underway. 

Results 

Various  tests  have  been  performed  to  determine  the  settings  that  produced  the  most 
precise  and  least  noisy  attenuation  map.  Initial  tests  varied  the  ultrasound  power,  gain  setting, 
and  whether  or  not  the  region  of  interest  should  be  two-dimensional  (space  averaged)  or  three- 
dimensional  (volume  averaged).  For  these  preliminaiy  experiments,  visual  comparison  of  the 
attenuation  images  versus  the  anticipated  attenuation  image  (obtained  from  examining  the  B- 
mode  image)  was  used  to  determine  the  optimal  settings. 

Preliminaiy  results  showed  no  benefit  from  the  3D  volume  averaging.  During  the  analysis 
of  the  attenuation  images,  it  was  noted  that  data  was  lost  in  areas  where  the  ultrasound  back 
scatter  was  too  low,  or  too  high.  For  any  single  gain  setting,  the  data  loss  was  significant. 
However,  based  on  a  visual  comparison,  other  gain  settings  appeared  to  produce  attenuation 
images  of  the  almost  the  same  quality,  but  with  data  losses  in  different  regions.  These 
observations  resulted  in  the  multiple  gain  setting  technique  discussed  above.  This  technique 
resulted  in  images  with  significantly  less  data  loss. 

Figure  1 1  presents  a  B-mode  image  of  a  breast  biopsy  phantom.  This  image  was  taken  at 
the  same  time  that  attenuation  data  was  being  acquired  for  the  same  image  space.  Because  the 
attenuation  slope  value  is  assigned  to  a  point  in  the  center  of  the  ROI,  and  the  because  of  the  size 
of  the  ROI,  0.5  cm  of  attenuation  data  was  lost  on  each  side  of  the  B-mode  image.  For  the  same 
reason,  0.85  cm  was  lost  from  the  top  and  bottom  of  the  image.  Figure  12  shows  the 
corresponding  attenuation  map  for  the  region  indicated  in  Figure  1 1 . 

Discussioii 

Attenuation  maps  are  not  a  new  application. 

Walach  and  others  have  produced  attenuation  maps  in  the 
past  [12,13,23].  However,  this  application  is  the  first 
three-dimensional  mapping  system  we  have  been  able  to 
find.  This  application  has  limitations  similar  to  the  other 
maps  that  have  been  produced.  One  such  limitation  is  the 
resolution  of  the  attenuation  map. 

Most  attenuation-calculation  techniques  require  a 
large  ROI  in  the  depth  direction.  The  literature  has 
recommended  that  the  log  spectral  difference  method  be 
used  for  global  attenuation  estimation  (greater  than  5.0  cm) 

[14,22].  However,  Kuc  has  theorized  that  this  method  can 
be  used  for  smaller  regions,  but  that  the  size  of  the  ROI  has 
a  lower  limit.  As  the  length  of  the  ROI  decreases  below 
this  limit,  the  noise  associated  with  the  application 
increases  dramatically  [26].  As  can  be  seen  from  Figure 
12,  the  resolution  is  poor  compared  to  the  B-scan  image. 

In  addition,  the  attenuation  image  only  partially  matches 
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Figure  12:  Attenuation  image.  This  image 
corresponds  to  the  highlighted  region 
in  Figure  2.  Values  of  attenuation  are 
inl/(cmdBMHz). 
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the  expected  image  (from  examining  the  B-mode  image  in  Figure  1 1). 

Inhomogeneous  regions  are  another  problem  with  this  technique.  Several  interfaces  and 
varying  structures  were  present  in  the  phantom  that  was  used  for  calculation.  The  development 
of  this  and  other  estimation  techniques  were  based  on  the  tissue  being  homogeneous  within  the 
ROI  [14,17,18,21,22].  This  problem  was  initially  ignored  in  this  study  in  order  to  discover  how 
the  method  would  handle  interfaces  and  different  tissues  within  the  ROI..  The  bright  circle  in 
Fig.  1 1  is  a  water  filled  sphere,  1.2  cm  diameter.  Note  diat  while  this  sphere  is  smaller  than  the 
ROI  it  does  show  up.  This  is  taken  as  a  sign  that  this  method  may  be  used  to  calculate  the 
attenuation  in  small  tissue  structures. 

Another  source  of  error  in  the  three-dimensional  map  is  the  location  error  of  the 
transducer.  In  tests,  the  largest  error  in  positioning  the  transducer  found  was  0.04  mm.  The 
beam  width  of  this  transducer  at  the  focal  point  is  1 .4  mm.  Therefore,  the  positioning  error  is  not 
considered  a  significant  source  of  error. 

Regardless  of  the  errors  involved,  the  attenuation  estimation  mediod  should  be  successful 
if  it  can  appropriately  predict  the  relationship  between  die  attenuation  and  the  various  structures 
within  the  phantom  or  anatomy.  (i.e.  which  regions  have  “higher”  attenuation  and  which  regions 
have  lower  attenuation.)  This  kind  of  an  approach  can  be  taken  because  of  the  iterative  nature  of 
the  inverse  model  that  will  use  these  values. 


SHbsystem  2:  Experimental  SAR  Measarcment 

•  Thermocouple  based  beam  plots  have  been  performed  to  calibrate  the  ultrasound 
power  field  in  water  for  specific  transducers. 

•  Experiments  have  been  conducted  in  water  and  will  begin  shortly  using  cast  agar 
phantoms. 

•  Preliminaiy  in-vivo  experiments  were  conducted  on  a  dog. 

•  A  finite  difference  model  has  been  constructed  to  help  interpret  the  experimental 
measurements 

To  accomplish  the  inverse  optimization  of  the  theoretical  model,  the  results  of  the  model 
must  be  compared  against  actual  values  that  are  determined  by  experimental  means.  The 
ultrasound  propagation  model  predicts  the  incident  power  at  all  points  in  the  modeled  region,  as 
well  as  the  absorbed  power.  At  a  given  point,  these  two  quantities  are  related  through  the  local 
attenuation.  Thus,  the  model  can  be  optimized  by  comparing  the  propagation  model  results  to 
either  the  incident  power  or  the  absort^d  power  at  several  points  in  the  model.  The  experimental 
values  needed  for  optimization  can  be  determined  using  the  same  thermocouples  which  are 
inserted  for  observation  of  steady  state  temperature  during  clinical  hyperthermia  treatment. 


Experimental  measurement  of  the  power  deposited  in  tissue  at  a  specific  point  is  not 
entirely  straight  forward.  Researchers  note  that  the  interaction  of  the  thermocouple  probe  with 
the  ultrasound  results  in  temperature  artifact.  [24]  This  artifact  results  fitim  two  sources:  viscous 
heating  at  the  interface  between  the  probe  surface  and  the  tissue,  and  absorption  of  acoustic 
energy  by  the  probe  its  self.  The  viscous  heating  is  a  surface  phenomena  and  results  in  a  small 
artifact  with  a  short  time  constant  causing  a  rapid  initial  temperature  rise  when  acoustic  power  is 
initiated.  This  artifact  is  usually  only  noticed  when  conducting  tests  using  bare  thermocouples. 
Probe  sheaths  of  steel,  fused  silica,  and  various  plastics  have  been  shown  to  absorb  ultrasound 
power  at  a  higher  rate  than  tissue.  Plastics,  the  most  common  sheath  materials  used  in 
hyperthermia  treatment,  absorb  ultrasound  energy  at  a  significantly  higher  rate  than  tissue 
resulting  in  a  significant  artifact  that  is  exhibited  over  a  longer  time  scale  than  the  viscous 
artifact.  The  temperature  rise  in  the  tissue,  a  much  larger  volume  than  the  sheath,  occurs  over  an 
even  longer  time  scale. 

The  combined  result  of  these  effects  is  shown  in  Figure  13.  In  Figure  13(a),  a  bare 
thermocouple  was  submersed  in  a  water  bath  and  insonified  for  0.6  seconds.  Note  the  rapid 
initial  jump  followed  by  a  relatively  linear  region.  The  initial  jump  is  caused  by  viscous  heating 
at  the  thermocouple/water  interface  and  the  linear  region  is  established  when  local  conduction 
between  the  thermocouple  and  surrounding  water  becomes  significant.  At  longer  times  the 
temperature  rise  will  become  non-linear  as  larger  scale  conduction  and  convection  become 
significant.  In  Figure  1 3(b)  the  experiment  was  repeated  with  a  sheathed  thermocouple.  Note 
that  the  temperature  rise  is  much  more  significant,  rising  roughly  20  degrees  in  0.6  seconds 
compared  to  less  than  one  degree  for  the  bare  thermocouple.  This  large  jump  is  attributed  to  the 
power  absorption  of  the  sheath.  Any  viscous  heating  that  occurs  is  clearly  lost  in  diis  rapid 
temperature  rise.  The  slope  of  this  rise  is  indicative  of  the  power  absorption  in  the  sheath. 
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If  the  power  is 
applied  for  a  longer  time, 
the  rapid  initial  rise  of  the 
sheathed  thermocouple  is 
balanced  by  conduction 
between  the  sheath  and  the 
surrounding  water.  This 
effect  is  shown  in  Figure 
14.  The  second  linear 
regime  achieved  represents 
the  temperature  rise  of  the 
media  surrounding  the 
thermocouple. 

In  diese  experiments 
the  thermocouples  were 
immersed  in  water,  a  very 
poor  absorber  of  ultrasound 
energy.  Thus,  essentially 
all  of  the  temperature  rise  of 
the  water  was  due  to 
conduction  from  the 

thermocouple  sheath.  In  tissue  we  expect  larger  absorption  directly  by  the  tissue.  An  actual 
experiment  conducted  in  tissue  will  produce  results  similar  to  those  shown  in  Figure  15.  Here  a 
sheathed  thermocouple  was  inserted  in  the  muscle  tissue  of  the  thigh  of  a  dog.  The  focus  of  the 
ultrasound  transducer  was  centered  on  the  thermocouple  and  the  area  was  insonified  for  0.6 
seconds  at  several  different  power  levels.  Note  that  in  each  case  the  power  was  much  lower  than 
that  used  in  the  experiment  from  Figure  14,  evidenced  by  the  relatively  low  temperature  rise.  In 
future  experiments  we  intend  to  calibrate  the  plastic  sheath  and  use  it  as  an  incident  power 
sensor,  allowing  the  ultrasound  power  incident  at  the  thermocouple  to  be  measured. 


Modeling  the  Thermocouple  Behavior 

In  the  absence  of  conduction,  power  is  related  to  heat  generation  by  the  following 
equation; 


dT  , 

~  ^  Equation  1 

where  p  is  the  density,  Cp  is  the  specific  heat,  a  is  the  absorption  coefficient,  and  I  is  the  local 
intensity  of  the  ultrasound  field. 


If  a  plastic  thermocouple  probe  and  surrounding  tissue  at  uniform  temperature  is  suddenly 
insonated  by  ultrasound  of  constant  intensity,  the  probe  and  tissue  will  begin  to  rise  in 
temperature.  The  plastic  probe  will  heat  more  rapidly  than  the  tissue  because  it  has  a  larger 
absorption  coefficient  than  the  tissue.  For  a  short  period  of  time  immediately  following 
insonation,  conduction  of  heat  will  be  small  and  the  rate  of  temperature  rise  is  dependent  solely 
on  the  absorption  coefficient  of  the  medium  being  insonated.  If  the  time  constant  of  the 
thermocouple  is  sufficiently  small  the  thermocouple,  being  surrounded  by  plastic,  will  reflect  the 


20 


41 

40.5 

40 

39.5 
^  41 

Si  40.5 
m  40 

3  39.5 

<  41 

CC 

^  40.5 

2  40 
UJ 

^  OQ  C 

- r—  1 — 

Impulse  test  of  dog  at  '35%''  power 
rate  of  heating  is  0.68  C/s 

1  2 

3  4  5  6 

Impulse  test  of  dog  at  '40%'  power 
rate  of  heating  is  0.77  C/s 

1  2 

3  4  5  6 

Impulse  test  of  dog  at  "45%"  power 
rate  of  heating  is  1.03  C/s 

^  dy.5 

41 

1  2 

- 1 - - — 1 - ; - ”  1 

3  4  5  6 

40.5 

XV 

Impulse  test  of  dog  at  '50%'  power 

40 

^  ^rate  of  heating  is  1 .39  C/s 

39.5 

>  .  «  -i- . 

TIME  (seconds) 

Figure  15:  Four  sequential  tests  of  a  plastic  sheathed  thennocouple  inserted  in  the  thigh  of  a 
_ dog.  _ 

temperature  of  the  plastic.  If  the  absorption  characteristics  of  the  thermocouple  probe  can  be 
sufficiently  characterized,  the  intensity  of  incident  ultrasound  can  be  determined  by  the 
temperature  vs.  time  response  of  the  ^ermocouple  to  a  suddenly  applied  ultrasound  field. 

The  radial  distribution  of  a  focussed  ultrasound  beam  can  be  closely  estimated  as  a 
Gaussian  function.  The  following  equation  describes  the  radial  intensity  distribution: 

A'*)  =  4iax^  ^  Equation  2 

where  Imax  is  the  intensity  at  r  =  0  and  p  is  a  measure  of  the  beam  width  as  determined  by  a  least 
squares  error  curve  fit.  [24] 

If  the  ultrasound  is  applied  for  a  short  duration  to  prevent  significant  conduction  effects 
then  the  temperature  distribution  produced  by  the  absorption  of  the  ultrasound  will  also  be 
Gaussian.  Using  the  theory  put  forth  by  Parker,  the  time  dependent  temperature  at  the  center  of 
the  focal  region  can  be  described  by  the  following  equation: 

^(0  =  Tnax  1)  Equation  3 

where  k  is  the  heat  conduction  coefficient  of  the  tissue  and  Jmax  is  the  maximum  temperature 
attained  at  the  center  of  the  focal  region  at  the  end  of  ultrasound  application. 
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At  the  tennination  of  the  momentaiy  application  of  ultrasound  the  plastic  sheathing  of  the 
thermocouple  probe  will  be  at  a  higher  temperature  than  die  surrounding  tissue.  However, 
because  the  heat  capacity  of  the  thermocouple  is  small  compared  to  the  heat  capacity  of  the 
tissue,  excess  heat  stored  in  the  probe  will  diffuse  into  the  tissue  and  the  probe  temperature  will 
reach  equilibrium  with  the  surrounding  tissue  and  follow  the  temperature  decay  of  the 
surrounding  tissue.  Curve  fitting  die  decay  of  the  tissue  temperature  by  the  least  squares  error 
method  the  decay  can  be  backtracked  and  the  value  of  Tnwx  determined.  Tmax  and  the  duration  of 
insonation  are  used  with  equation  (1)  and  the  intensity  determined  by  the  rate  of  heating  to 
determine  the  absorption  coefficient  of  the  tissue.[24] 

It  has  been  determined  by  previous  investigators  that  absorption  accounts  for  90-100%  of 
attenuation  in  soft  tissue.[25]  Therefore  it  is  not  completely  accurate  to  compare  attenuation 
determined  by  B-mode  imager  and  absorption  determined  by  thermocouple  probe  but  an  error 
less  than  10%  is  not  unreasonable  when  compared  to  other  errors  likely  to  exist. 

A  program  was  created  on  the  hospital  treatment  system  that  positions  the  transducer, 
delivers  power  for  a  period  of  time,  and  records  the  temporal  temperature  data  fi'om  the 
thermocouple  probes.  This  program  is  able  to  attain  temperature  data  fi'om  several  seven 
junction  probes  every  3/100  sec.  To  test  the  program  a  10  cm,  2MHz  transducer  was  mounted  to 
the  positioning  gantry  and  an  array  of  thermocouples  was  placed  in  a  water  bath  above  the 
treatoent  system.  The  thermocouples  were  located  with  the  existing  thermocouple  locate 
routine.  It  was  found  that  the  positioning  system  was  only  able  to  locate  the  thermocouples  with 
i  2-4cm  accuracy.  A  radial  beam  intensity  plot  was  attempted  begiiming  1  cm  away  fi’om  the 
center  of  focus  and  traversing  through  the  focus  at  2mm  increments.  Because  of  the  tight  beam 
focus  and  lack  of  positioning  accuracy  it  was  found  that  there  was  very  little  correlation  between 
expected  position  and  beam  intensity. 

To  characterize  the  beam  shapes  of  the  available  transducers,  an  array  of  epoxy  covered 
thermocouples  was  created  and  oriented  perpendicular  to  the  beam  axis.  The  array  was  moved 
through  the  beam  transverse  to  the  beam  axis  in  Inun  increments  and  the  steady  state  temperature 
was  recorded  for  the  thermocouples  at  several  different  distances  fi-om  the  face  of  the  transducer. 
This  procedure  was  followed  for  the  focussed  10cm  diameter  2MHz,  focussed  7cm  diameter 
500kHz,  unfocussed  square  4cm  IMHz,  and  unfocused  3cm  diameter  4MHz  transducers.  The 
unfocussed  transducers  were  not  characterizable  because  of  the  irregular  intensity  consistent  with 
the  long  near  field  of  the  beams.  The  1 0cm  and  7cm  diameter  transducers  had  a  Gaussian  profile 
in  the  radial  direction  at  the  focus.  Least  squares  error  curve  fit  produced  a  Gaussian  parameter 
of  0.029  cm'  for  the  10  cm  diameter  transducer  and  a  Gaussian  parameter  of  .42  cm'^  for  the 
7cm  diameter  transducer. 

The  0.029  cm  ^  Gaussian  parameter  at  the  focus  of  the  10cm  diameter  transducer  creates  a 
beam  width  of  around  6mm.  which  is  too  narrow  for  the  accuracy  of  the  hospital  positioning 
system  and  die  steep  temperature  gradients  created  by  the  narrow  focus  decrease  the  temporal 
window  during  which  conduction  is  insignificant.  The  0.42  cm’^  Gaussian  parameter  at  the  focus 
of  the  7  cm.  diameter  transducer  creates  a  beam  width  larger  than  2  cm.  This  larger  beam  width 
will  create  a  less  steep  temperature  gradient  and  lengthen  the  temporal  window  for  the  rate  of 
heating  measurement.  The  longer  length  of  focus  of  the  7  cm.  transducer  is  also  beneficial  in 
avoiding  error  in  the  Z-direction. 
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Finite  difference  code  was  written  to  model  a  thermocouple  probe  embedded  in  tissue. 
The  finite  difference  model  will  be  used  to  verify  the  assumptions  and  clarify  the  unknowns  of 
the  experimental  technique.  Verification  of  the  1-D  model  is  in  process  after  which  a  2-D  and/or 
3-D  model  will  be  created. 

A  considerable  amount  of  phantom  testing  with  the  7cm.  diameter  transducer  will  begin 
shortly.  The  phantom  tests  will  be  designed  to  prove  the  accuracy  of  the  rate  of  heating  and 
thermal  pulse  decay  technique  in  combination  with  the  probe  geomeny  used  in  clinical 
hyperthermia  treatment.  After  verification  the  procedure  will  be  integrated  with  the  other 
subsystems  of  die  inverse  estimation  technique  to  provide  geometric  data  for  the  thermocouple 
probes  and  local  measured  power  and  absorption  values. 


Subsystem  3:  Model  Development 

•  The  parabolic  model  has  been  selected  for  use  in  regions  of  low  acoustic  contrast.  In 
these  regions  the  model  is  quite  fast  and  loses  little  accuracy  when  compared  to  the 
Green’s  Function  Model.  The  Green’s  Function  Model  can  be  used  in  high  contrast 
areas  (close  to  ribs,  etc.) 

•  An  interface  has  been  built  for  the  model  which  allows  the  user  to  specify  the  volume 
modeled,  the  acoustic  power  input  as  a  2D  complex  pressure  function,  and  the  model 
parameters;  density,  sound  velocity,  and  absorption  at  each  node  within  the  volume. 


Based  on  the  results  of  the  2nd  year,  the  parabolic  model  was  selected  for  use  in  regions 
of  low  acoustic  contrast.  In  these  regions  the  model  is  quite  fast  and  loses  little  accuracy  when 
compared  to  the  Green’s  Function  Model.  If  necessary,  the  Green’s  Fimction  Model  can  be  used 
in  high  contrast  areas  (close  to  ribs,  etc.)  The  Green’s  Function  Model  could  be  used  everywhere 
resulting  in  a  significantly  slower,  slightly  more  accurate  result. 


During  the  third,  year  work  on  the  model  focused  on  the  development  of  a  practical  and 
convenient  user  interface  that  allows  the  model  to  be  adapted  to 
different  transducers  and  tissue  geometries.  Toward  this 
purpose,  an  interface  was  built  for  the  model  which  allows  the 
user  to  specify  the  volume  modeled,  the  acoustic  power  input  as 
a  2D  complex  pressure  function,  and  the  model  parameters, 
density,  sound  velocity,  and  absorption  at  each  node  within  the 
volume.  The  model  interface,  as  described  below  is  complete. 

The  interface  allows  the  user  to  specify  the  dimensions, 

X,  Y,  Z  of  a  rectangular  volume  that  will  be  modeled.  The  size 
of  the  volume  to  be  modeled  is  limited  by  the  requirement  for 
2  nodes  per  wavelength  at  the  transducer  frequency  under 

investigation  and  by  the  maximum  number  of  nodes  that  the  model  can  handle.  At  the  current 
time,  the  model  can  handle  1,000,000  nodes.  Thus,  a  volume  of  100  by  100  by  100  nodes  could 
be  modeled.  At  1 .5  MHz,  this  would  be  a  5  cm  cube.  It  should  be  noted  that  larger  volumes  can 
be  modeled  by  slicing  the  volume  into  1,000,000  node  rectangles  and  modeling  each 
sequentially. 


Within  the  modeled  volume 
the  tissue  is  characterized  by 
specifying  the  density,  sound 
velocity,  and  acoustic  absorption  at 
each  node.  As  a  result  the  interface 
is  very  general.  There  is  much  more 
flexibility  in  this  interface  than  is 
currently  required  by  the  inverse 
optimization  scheme.  If  we  later 
desire  to  optimize  the  sound  velocity 
field  within  the  tissue,  or  the  density 
field,  no  modifications  to  the  model 
will  be  required. 

Transducer  Model 

If  the  transducer  were 
modeled  as  part  of  the  above 
described  volume,  it  would  be 
difficult  to  fit  the  desired  tissue 
volume  into  the  limited  model  volume.  This  problem  is  addressed  by  solving  for  the  complex 
pressure  field  created  by  the  transducer  being  modeled  and  applying  this  field  at  one  surface  of 
the  model  volume. 

The  3D  complex  pressure  field  is  first  found  for  the  transducer  in  water.  This  is  done  by 
matching  a  numerical  solution  to  the  results  of  calibration  experiments  where  the  temperature  rise 
of  epoxy  coated  thermocouples  is  measured  in  a  water  bath.  Using  this  data,  as  long  as  the  edge 
of  the  modeled  volume  extends  beyond  the  skin  into  the  water  bath,  the  effects  of  the  transducer 
can  be  modeled  accurately  by  finding  the  intersection  between  the  planar  surface  of  the  model 
and  the  3D  pressure  field.  Note  that  a  single  calibration  allows  the  transducer  to  be  positioned 
anywhere  with  respect  to  the  model  volume. 

With  this  user  interface  in  place,  the  model  can  be  treated  as  a  black  box,  where  the 
optimization  scheme  is  ignorant  of  the  details  of  the  model  within  the  box,  but  simply  specifies 
inputs  and  collects  the  output  from  the  model.  Once  the  transducer  has  been  calibrated  (or 
several  transducers),  the  user  simply  specifies  the  transducer  location  and  orientation,  and  the 
nodal  acoustic  properties  that  were  identified  in  Subsystem  1 .  The  model  then  outputs  the 
acoustic  power  deposited  at  each  node  within  the  model  volume. 

Subsystem  4:  Inverse  Technique 

•  Basic  algorithm  development  is  complete.  Optimization  runs  have  been  completed  on 
simple  geometries. 


Fig  1.  surbce  plot  of  function  space  (-.05,  -.03) 
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Figure  16:  This  is  a  surface  plot  of  an  example  error  function 
space.  The  optimization  scheme  is  searching  for  the 
minimum  in  this  function. 


Inverse  Method  for  Optimizing  Tissue  Attenuation  Values 
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In  the  original  proposal,  it  was  suggested  that  three  parameters  could  be  allowed  to  vaiy 
between  tissue  types.  The  numerical  model  developed  has  the  capability  to  change  all  three 
parameters.  However,  in  initial  trials,  it  is  assumed  that  speed  of  sound,  (c)  and  impedance,  (Z) 
are  known  constants  for  the  tissues  with  which  we  will  be  dealing.  Consequently,  attenuation  is 
the  only  parameter  that  changes  between  tissue  types  and,  thus,  will  be  the  only  parameter 
optimized. 

The  method  discussed  in  the  proposal  uses  the  Jacobean  matrix  to  adjust  the  initial 
parameter  values  at  model  nodes  and  iterates  until  the  error  function  is  minimized.  (Here  error  is 
some  weighted  sum  of  the  differences  between  model  predictions  and  experimental 
measurements.)  This  method  was  used  in  initial  trials  and  then  compared  to  other  methods.  A 
second  method,  the  Conjugate  Gradient  method,  was  also  evaluated. 

Trials  were  completed  using  each  of  these  methods  to  search  for  minimum  error  points  in 
the  error  function  for  simple  function  spaces.  Figure  16  shows  an  example  where  only  two 
parameters  are  varied.  Given  an  initial  value  for  the  two  parameters,  each  method  iteratively 
searches  for  the  optimum  (minimum)  point  in  the  function.  Essentially  the  only  difference 
between  the  various  optimization  schemes  is  the  efficiency  with  which  they  “slide  down  the  hill” 
to  find  the  minimum  value.  This  efficiency  is  based  on  the  number  of  function  calls  (forward 
model  iterations)  required  to  produce  one  optimization  iteration,  and  on  the  number  of 
optimization  iterations  required  to  reach  the  minimum  error.  Figure  17  shows  trials  using  both 
the  Jacobian  and  the  Conjugate  Gradient  methods  for  optimization. 


Number  of  function  calls  (one  ooint  oer  iteration) 


I 


Number  of  (unction  calls  (one  point  per  iteration) 


Figure  17:  This  shows  the  rate  of  convergence  for  several  versions  of  the  Conjugate  Gradient  (left)  and  the 
Jacobian  (right)  based  optimization  schemes. 
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Fourth  Year  Accomplishments: 

Subsystem  1;  Geometry  AcauisitioB  and  Atteuuatioa  Measurement 

•  Verification  tests  for  clinical  system  for  3D  attenuation  mapping. 

The  hardware  and  software  for  the  attenuation  mapping  system  were  completed  by  the 
end  of  the  third  year  of  the  project.  During  the  fourth  year,  validation  tests  were  performed. 
These  results  will  be  reported  quite  briefly  here.  Detailed  results  can  be  found  in  [27]. 

In  this  woric  a  modified  ultrasonic  B-mode  imaging  system  was  combined  wiA  a  two-axis 
motion  device,  making  it  possible  to  collect  a  series  of  raw,  reflected  ultrasonic  wave  forms  that 
represent  a  three-dimensional  volume.  The  attenuation  coefficient  for  a  specific  region  of 
interest  (ROI)  was  estimated  using  the  log  spectral  difference  technique.  In  order  to  reduce  the 
level  of  noise  in  the  attenuation  estimate,  a  region  of  interest  was  specified  within  which  the 
estimates  were  averaged.  Using  experimental  data  collected  fiorn  graphite/agar  phantoms  the 
minimum  size  ROI  that  will  produce  a  reliable  signal  has  been  investigated.  The  results  show 
that  a  2D  ROI,  where  information  is  taken  from  a  region  within  a  B-scan,  must  be  approximately 
2.5  cm  in  the  depth  direction  by  1  cm  in  width  to  produce  a  reasonable  signal.  Information 
combined  from  adjacent  B-scans  into  a  3D  ROI  can  reduce  the  apparent  size  of  this  region  within 
the  2D  plane  by  improving  the  lateral  resolution.  In  this  case,  the  minimum  ROI  was  identified 
as  2.5  cm  (depth)  by  5  mm  by  5  mm. 

Methods 

A  clinical  B-mode  imaging  system  and  a  two-axis  motion  machine  were  used  to  acquire 
reflected  ultrasoimd  signals  from  a  graphite/agar  phantom  in  three  dimensions.  Then,  the  log 
spectral  difference  technique  was  used  to  calculate  the  local  attenuation  coefficient.  The 
accuracy  of  the  attenuation  measurements  was  verified  with  a  force-balance  calculation  of 
attenuation  made  from  the  phantom.  Details  of  the  system  and  the  attenuation  calculation 
techniques  have  been  presented  in  previous  reports  and  will  not  be  presented  here.  The 
following  sections  present  the  methods  used  to  acquire  and  analyze  the  results. 

Data  Acqubition 

Preprocessed,  reflected  ultrasound  signals  were  acquired  from  a  modified  clinical 
imaging  system  and  a  two-axis  motion  system.  In  order  to  produce  3D  results,  multiple  slices  of 
a-lines  were  acquired  on  a  symmetric  2D  grid.  Each  slice  of  data  corresponded  to  the  B-mode 
image  produced  by  that  data.  Individual  a-lines  acquired  within  each  slice  were  spaced  1  mm 
apart  (x  direction).  Slices  were  spaced  1  mm  apart  in  the  direction  orthogonal  to  the  spacing 
between  a-lines  within  a  slice  (y  direction).  The  final  dimension  came  from  the  reflected  signal 
itself  (z  direction  or  depth).  Acquisition  of  the  signals  in  this  manner  was  accomplished  by 
fixing  the  imaging  transducer  to  a  two-axis  motion  system,  moving  the  transducer,  and  acquiring 
the  signals. 

The  signals  were  acquired  while  the  transducer  was  scanned  across  a 
homogeneous  graphite/agar  phantom.  A  layer  of  water  was  used  as  the  acoustic  pathway  to 
couple  the  imaging  transducer  and  the  phantom.  In  order  to  produce  a  standard  to  compare  with 
the  attenuation  measurement  system  results,  the  same  graphite/agar  solution  used  to  make  the 
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imaging  phantom  was  also  used  to  make  two  force  balance  phantoms.  Each  force  balance 
phantom  had  a  different  constant  thickness.  Therefore,  an  attenuation  coefficient  could  be 
calculated  from  force  balance  tests.  In  addition,  the  calculated  coefficient  from  the  force  balance 
could  be  used  to  determine  the  accuracy  of  the  attenuation  measurement  system. 

Data  Analysis 

The  log  spectral  difference  method  was  used  to  estimate  the  slope  of  the  attenuation  to 
frequency  relationship  [14, 15,  22, 4].  This  method  utilizes  the  ratio  of  the  power  spectrum  from 
a  shallow  and  deep  segment  of  a  reflected  ultrasound  signal  to  compute  the  attenuation  slope  for 
that  region.  Once  the  slope  has  been  calculated,  the  normalized  attenuation  coefficient  can  be 
computed  by  using  Equation  4.  In  this  equation  is  the  normalized  attenuation  coefficient 


P= 


6) 

2D 


Equation  4 


(normalized  to  depth  and  frequency),  a)  is  the  slope  of  the  log  spectral  difference,  and  D  is  the 
distance  between  the  shallow  and  deep  segments  and  is  equal  to  two-thirds  of  the  specified  ROI 
depth.  In  order  to  compute  an  attenuation  coefficient,  the  relationship  between  attenuation  and 
frequency  is  assumed  to  be  linear,  and  at  zero  frequency  the  coefficient  is  assumed  to  be  zero. 
Therefore,  the  calculation  of  the  attenuation  coefficient  (o)  can  be  easily  calculated  by  using 
Equation  5,  where /is  the  frequency.  Since  the  attenuation  coefficient  is  frequency  dependent, 
all  values  reported  in  this  paper  are  at  a  frequency  of  1.0  MHz.  This  method  was  used  in 
combination  with  signal  averaging  and  a  moveable  ROI  to  produce  an  attenuation  map  or  image. 

To  produce  the  attenuation  map,  a  ROI  of  a  determined  size  was  moved  throughout  the 

^  Equation  5 

three-dimensional  data  space  until  an  attenuation  map  was  produced.  The  dimensions  of  the  ROI 
were  varied  and  tested.  Three  different  ROI  settings  were  used  to  test  the  attenuation  calculation 
ability  of  the  system.  The  first  setting  was  a  2D  ROI  (Area  map)  where  the  width  (x  direction)  of 
the  ROI  was  1  cm.  This  produced  an  attenuation  calculation  that  was  averaged  over  10  a-lines. 
The  second  setting  was  a  small  3D  ROI  (Vol3  map)  where  the  width  was  3  mm  and  the  height  (y 
direction)  was  3  mm.  This  setting  produced  a  calculation  that  was  averaged  over  nine  a-lines. 
The  final  setting  was  a  3D  ROI  (Vol5  map)  that  was  slightly  larger.  The  x  and  y  dimensions  of 
this  setting  were  5  mm  by  5  mm.  This  setting  produced  a  calculation  that  was  averaged  over  25 
a-lines.  At  each  of  the  settings,  ROI  depths  (z  direction)  of  1.0,  1.5,  2.0,  2.5,  3.0,  and  3.5  cm 
were  tested. 

Once  the  ROI  size  was  selected,  it  was  moved  throughout  the  data  space  incrementally  by 
1 .0  mm  in  the  x  and  z  directions.  The  attenuation-slope  value  fiir  each  calculation  was  assigned 
to  a  point  in  the  center  of  the  ROI,  thus  producing  an  attenuation  map.  In  order  to  relate  die 
results  of  the  three  settings,  the  ROI  position  was  set  such  that  the  map  produced  from  each  ROI 
setting  represented  the  same  2D  space.  Therefore,  Vol5  map  required  two  data  slices  on  either 
side  of  the  center  slice,  Vol3  map  required  one  slice  on  either  side,  and  Area  map  did  not  require 
any  additional  slices. 
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Force  Balance  Comparison 


In  order  to  validate  the  attenuation  measurements  made  by  the  attenuation  mapping  system,  a 
force  balance  station  was  assembled.  This  station  used  a  reflecting  target,  an  ultrasound 
transducer  (different  from  the  imaging  transducer),  and  a  scale  to  measure  the  attenuation.  The 
reflecting  target  was  used  as  a  transducer  to  convert  acoustic  energy  into  a  mechanical  force, 
which  was  read  by  the  scale.  Two  slabs  of  the  graphite/agar  phantom  tested  in  the  attenuation 
mapping  system  were  tested  in  the  force  balance.  Each  slab  had  a  different  thickness  (</i  and  di). 
The  force  or  weight  ratio  (wi/wa),  as  read  on  the  scale,  is  the  same  as  the  power  ratio  of  the 
ultrasound  wave  as  it  passes  through  the  different  distance  of  the  same  material.  Therefore,  the 
weight  ratio  was  used  to  compute  the  attenuation  coefficient  at  the  frequency  of  the  transducer  by 
using  Equation  6.  The  attenuation  coefficient  calculated  from  the  force  balance  was  then  used  to 


-ln(-^) 

calculate  the  attenuation  coefficient  at  1.0  MHz  and  used 
mapping  system. 


Equation  6 

as  a  standard  for  the  attenuation 


Rcsalts 


This  study  is  focused  on  determining  the  accuracy  and  repeatability  of  attenuation 
measurements  when  the  ROI  dimensions  are  changed.  Three  sets  of  were  acquired  and 
analyzed  to  produce  an  attenuation  map  for  each  ROI  setting  and  depth.  From  this  map,  an  area 
of  interest  was  selected  where  the  entire  ROI  for  the  calculations  was  within  the  phantom  (there 
were  no  interfaces  within  the  ROI).  The  area  of  interest  was  approximately  3  cm  wide  and  1  cm 
deep  and  represented  approximately  300  computed  values  of  attenuation.  For  each  test,  the 
attenuation  coefficients  within  the  area  of  interest  were  averaged.  The  mean  values  were  then 
used  to  determine  the  accuracy  and  repeatability  of  the  measurement  system. 

The  results  show  that  the  accuracy  of  the  attenuation  measurement  system  is  most 
dependent  on  the  depth  of  the  ROI  (see  Figure  18).  For  all  test  settings,  an  ROI  deptfi  greater 
than  or  equal  to  2.5  cm  produced  accurate  estimations  of  the  attenuation  coefficient  for  the 
graphite/agar  phantom.  However,  histogram  plots  of  the  attenuation  coefficient  from  an  area  of 
interest  show  that  the  Vol5  map  setting  (5mm  by  5mm  in  the  x  and  y  directions)  is  the  best 
setting  for  calculating  an  attenuation  coefficient  (see  Figures  19-21).  In  addition  to  providing 
more  accurate  attenuation  calculations,  the  Vol5  map  setting  provides  improved  lateral  resolution 
within  a  2D-attenuation  image.  Therefore,  the  Vol5  setting  with  an  r6i  depth  of  2.5  cm  is  the 

best  setting  for  determining  the  attenuation  coefficient  of  small  regions  and  for  producing  an 
attenuation  map. 

Detailed  results  can  be  found  in  [27]. 


Standard  deviation  as  a  percentage 
of  the  force  balance  estimate 


Figure  18:  This  shows  the  decrease  in  the  standard  deviation  of  the  attenuation  predictions  as 
the  ROI  depth  is  increased. 
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Figure  19:  This  shows  the  changing  spread  of  the  histogram  of  attenuation  values  changing  as 
the  ROI  is  extended  in  the  depth  direction. 
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Figure  20:  This  shows  the  changing  spread  of  the  histogram  of  attenuation  values  changing 
the  ROI  is  extended  in  the  depth  direction. 
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•  Software  and  hardware  for  clinical  patient  geometry  acquisition  completed. 


During  the  fourth  year  the 
patient  geometry  acquisition  system 
was  completed.  This  required  the 
completion  of  the  Geometry 
Acquisition  and  Display  Soltware 
(GADS)  and  completion  of  the  patient 
registration  system. 

Figure  22  shows  the  clinical 
system  and  Figure  23  shows  the  basic 
hardware.  An  imaging  transducer  is 
mounted  on  the  gantry  that  moves  the 
treatment  transducer,  and  the  two 
transducers  are  moved  together. 


Figure  22:  Clinical  Hyperthermia  System. 
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Figure  23:  Schematic  of  completed  system. 
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As  the  motion  controller  moves  the  transducer  the  B-scan  imager  produces  images  that  are 
captured  by  the  Frame  Grabber. 

At  this  point  the  GADS  software  takes  over.  Image  sets  are  recorded  and  organized  along 
with  pertinent  patient  data.  In  GADS,  the  images  are  automatically  composited  using  the  known 
X,  y,  z  information  from  the  motion  controller.  In  this  example,  the  transducer  has  been  translated 
to  three  positions  that  lie  side  by  side,  in  the  plane  of  the  B-scan  image.  The  detailed  operation  of 
the  system  has  been  described  previously.  The  following  extended  Figure  (Fig.  24)  includes  a 
series  of  screen  captures  that  will  clarify  the  operation  of  the  software. 
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These  surfaces  separate  the  image  into  regions.  Each  region  is  assigned  a  tissue 
type,  which  assigns  model  parameters:  density,  speed  of  sound,  and  attenuation. 
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The  node  size  appropriate  to  the  model  is  selected  and  the  nodal  values  are  assigned. 
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The  results  are  displayed  visually  showing  the  effects  of  the  node  size. 


The  various  tissue  types  and  their  assigned  properties  are  recorded  in  a  table. 


Figure  24;  Operation  of  the  Geometiy  Acquisition  and  Display  Software  (GADS). 


Patient  Registration 


The  final  ^stem  developed  to  reliably  and  repeatably  position  the  patient  for  treatment 
uses  three  systems: 

1 .  Laser  Alignment 

2.  Thermoplastic  Immobilization 

3.  Ultrasonic  Maricers 

The  laser  alignment  system  uses  tattoos  on  the  patient's  skin  that  are  aligned  to  lasers 
permanently  positioned  above  and  to  the  side  of  the  treatment  table  (see  Figure  25).  Once 
positioned  using  the  laser  system,  the  patient  is  locked  into  place  using  a  thermoplastic  sheet. 
This  sheet  is  molded  to  die  patient's  botfy  during  the  first  visit  and  helps  return  the  patient  to  an 
identical  position  on  subsequent  visits  (see  Figure  26).  1/4  inch  diameter  disks  of  this  steel  shim 
material,  acting  as  ultrasonic  maikers,  are  placed  on  the  patient's  skin  near  the  treatment  region 
and  are  used  to  verify  the  location  of  soft  tissue  in  the  treatment  region.  Figure  27  shows  the 
appearance  of  the  markers  in  a  B>scan  image. 


Figure  26:  This  shows  the  design  schematic  for  the  Thermoplastic  Positioning  System 
along  with  the  completed  system. 


Figure  27:  This  shows  the  appearance  of  the  ultrasound  maricers  in  a  B-scan  image. 

There  are  two  possible  artifacts  associated  with  these  markers.  On  the  left, 
the  markers  create  shadows  where  little  ultrasound  energy  passes  through 
the  marker  to  tissue  below.  On  the  right,  the  marker  happens  to  be  normal 
to  the  transducer  face  and  multiple  reflections  fiom  the  marker  are  seen. 


•  Verification  tests  completed 

The  following  experiment  was  performed  in  order  to  estimate  the  accuracy  and 
repeatability  of  the  ultrasound  geometry  acquisition  system. 

Experimental  Setup 

Nine  five-millimeter  diameter  ultrasound-reflecting  markers  were  punched  out  of  a  sheet 
of 0.003-inch  thick  steel  shim  using  a  circular  paper  hole-punch.  The  circular  markers  were 
adhered  to  a  Gammex™  RMI-429  Ultrasound  Biopsy  Phantom  using  a  thin  coat  of  contact 
cement.  The  relative  positions  of  the  center  of  each  marker  were  measured  with  respect  to  the 
central  marker  (Marker  4)  using  the  three-axis  positioning  system  of  a  milling  machine.  The 
phantom  was  fixed  to  the  milling  machine  table  so  that  the  markers  were  on  the  superior  side  and 
a  pointed  mill  bit  was  visually  aligned  with  the  center  of  each  marker.  Figure  28a  shows  a 
schematic  top  view  of  the  markers’  relative  positions  on  the  biopsy  phantom  and  Figure  28b  plots 
the  measured  marker  positions  using  Marker  4  as  the  origin  of  the  coordinate  system.  Figure  29 
is  a  general  schematic  of  the  phantom  installed  in  the  test  apparatus. 

The  central  marker  (Marker  4)  lay  on  an  approximately  horizontal  surface  at  the  peak  of 
the  phantom.  Due  to  the  dome  shape  of  the  biopsy  phantom,  the  four  markers  closest  to  this 
central  marker  (Maricers  2, 3, 5,  and  6)  were  angled  approximately  30  degrees  off  the  horizontal 


plane.  Similarly,  the  markers  furthest  from  the  center  (Markers  1,  7,  8,  and  9)  were  tilted 
approximately  45  degrees  off  the  horizontal  plane. 

Results 

The  results  showed  that  reflectors  could  be  used  to  detect  patient  movements  with  an 
accuracy  of  approximately  1.5  millimeters.  The  repeatability  found  between  three  tests  was 
exceptional,  with  the  error  calculated  to  be  less  than  approximately  one  millimeter. 

Results  from  this  experiment  suggest  that  an  averaging  algorithm  that  would 
automatically  calculate  the  central  location  of  a  series  of  points  selected  on  a  particular  marker 
could  substantially  improve  the  accuracy  and  efficiency  of  the  technique.  The  information  could 
also  be  used  to  translate  the  treatment  plan  coordinates  to  the  new  patient  position,  saving 
operators  valuable  time  that  would  otherwise  be  spent  repositioning  the  patient.  A  similar  system 
could  also  be  used  to  aid  in  locating  thermocouple  catheters  prior  to  treatments  in  order  to 
expedite  the  process  of  finding  thermocouple  junctions  and  to  help  reduce  pain  often  associated 
with  this  search.  More  detailed  information  is  available  in  [28]. 


Figure  28:  a.)  Top  view  of  ultrasound-reflecting  maricers  on  the 
surface  of  the  biopsy  phantom,  b.)  Measured  maiker 
positions  with  respect  to  Marker  4. 
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Figure  29:  Schematic  of  linear  motion  apparatus,  phantom,  and  ultrasound 
transducer  in  the  water  tank. 


4“^ 


Subsystem  2;  Experimental  SAR  Measurement 

•  Probe  orientation  effects  on  the  SAR  measurements  were  investigated 
experimentally. 

A  re-examination  of  experimental  results  from  Year  3  and  a  series  of  new  tests  using  agar 
phantoms  showed  unexplained  variations  in  the  thermocouple  based  SAR  measurements. 
Continued  experiments  in  Year  4  document  and  explain  these  variations.  A  new  thermocouple 
probe  was  designed  to  minimize  these  variations. 

Figure  30  is  the  result  of  an  experiment  where  a  clinical  thermocouple  probe,  described 
below,  was  immersed  in  a  water  bath  and  each  of  the  seven  thermocouple  junctions  was  in  turn 
heated  by  placing  it  in  the  focus  of  an  ultrasonic  treatment  transducer  and  turning  the  transducer 
on  briefly.  While  the  ultrasound  power  and  relative  positioning  between  thermocouple  and 
transducer  were  held  constant,  the  rate  of  heating  measurement  varied  dramatically  from 
thermocouple  to  thermocouple  and  from  measurement  to  measurement  when  repeating 
experiments  on  the  same  thermocouple. 

In  order  to  explain  this  variation  (and  perhaps  eliminate  it)  a  series  of  experiments  were 
performed  where  the  angle  of  ultrasound  incidence  on  the  probe  was  varies  as  well  as  the  axial 
orientation  of  the  probe  (i.e.  the  probe  was  rotated  about  it's  axis).  This  experiment  and  results 
are  described  below. 
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LABORATORY  TESTING  APPARATUS 
Ultrasound  Power  Field  Generation 

A  diagram  of  the  laboratoiy  testing  apparatus  can  be  seen  in  Figure  3 1 .  The  input  power 
signal  was  generated  by  a  function  generator  and  amplified  by  an  rf  amplifier.  After 
amplification  the  power  signal  was  routed  through  a  power  pick  off  box  with  attached  power 
meter  to  measure  forward  and  reverse  electrical  power.  The  power  signal  continued  out  of  the 
power  pick  off  box,  passed  through  an  impedance  matching  box  and  terminated  at  the  ultrasound 
transducer.  Any  one  of  several  different  custom-built  ultrasound  transducers  can  be  chosen  to 
deliver  an  intensity  field  in  the  test  tank.  Table  2  is  a  list  of  the  transducers  available  for  use 


Table  2:  Available  Ultrasound  Transducers 

Shape  Diameter  or  Width  Foca 

- - - _ 

square  4 

square  4 

circular  3 

circular  7 

circular  10 

circular  10 


Focal  Length  (cm)  Frequency  (MHz) 


3-Axis  Positioner 


Rotating  Thermocouple 
Bradret 


Thermocouple  Probe 


Figure  3 1 ;  Diagram  of  experimental  system. 


Temperature  Measurement 


Thermocouple  output  was  measured  by  a  computer  connected  to  a  custom-built 
thermocouple  interface  box.  The  thermocouple  interface  box  can  rapidly  sample  as  many  as  16, 
7-junction  thermocouple  probes.  The  thermocouple  probe,  which  was  described  by  Anhault  and 
Hynynen,1992,  and  illustrated  in  Figure  32,  consists  of  seven  Manganin  Constantan 
thermocouple  junctions  spaced  at  one  centimeter  intervals  and  encapsulated  in  an  air  filled 
polyethylene  sheath  (o.d.  0.96-mm  and  i.d.  0.58-mm).  The  thermocouple  junctions  were 
manufactured  by  soldering  together  a  50-pm  Manganin  wire  and  50-pm  Constantan  wire.  The 
diameter  of  the  soldered  sensor  varies  but  is  approximately  0.3-mm  in  diameter.  A  variable 
temperature  isothermal  water  bath  was  used  to  calibrate  each  probe  using  a  two-point  calibration 
technique.  Calibration  data  for  each  thermocouple  was  stored  in  the  computer  and  used  during 
temperature  measurement. 

Mounting  and  Positioning  System 

The  three-axis  positioning  system  shown  in  Figure  29  was  placed  adjacent  to  the  water 
tank  so  that  the  positioning  arm  overiiung  the  tank.  The  L-shaped  transducer  mounting  bracket 
was  attached  to  the  positioning  arm  and  extended  towards  and  across  the  bottom  of  the  tank.  The 
transducer  was  fixed  to  the  transducer  mounting  bracket  and  leveled  so  that  the  ultrasound  field 
was  generated  towards  the  surface  of  the  water.  The  tank  was  filled  with  degassed  water  with 
oxygen  content  of  less  than  2-ppm.  During  early  ultrasound  experiments  bristle  brushes  were 
used  to  absorb  and  disperse  the  ultrasound,  however  it  was  apparent  from  test  results  that  some 
reflection  was  returning  from  the  surface  back  to  the  test  specimen.  A  reflecting  target  made  of 
closed  cell  foam  was  substituted  for  the  bristle  brushes  and  towels  were  draped  around  the  inside 
of  the  water  tank  to  completely  direct  and  absorb  the  excessive  ultrasound  intensity  away  from 
the  test  specimen.  Platforms  for  holding  phantoms  and  brackets  for  holding  thermocouples  in  the 
energy  field  were  fixed  to  a  steel  fi^e  that  encased  the  water  tank.  Several  different  apparatus 
were  used  during  various  tests  to  hold  and  position  the  test  specimens,  but  the  fixture  below  was 
used  in  the  experiments  described  here. 
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Figure  33:  This  is  the  thermocouple  orientation  fixture.  Using  this  fixture  it  is  possible  to 
rotate  the  axis  of  the  thermocouple  probe  with  respect  to  the  incident  ultrasound 
field,  or  to  rotate  the  probe  about  its  axis. 


Results 

When  the  axis  of  the 
probe  was  rotated  with  respect  to 
the  incident  power  field,  the 
result  shown  in  Figure  34  was 
produced.  This  figure  represents 
three  trials  as  the  probe  was 
rotated  to  a  position,  then  heated, 
then  rotated  again.  While  there  is 
some  scatter  between  repetitions, 
it  is  clear  that  the  measured 
power  drops  by  over  30%  once 
the  angle  of  incidence  reaches 
30°  away  from  normal  incidence. 
While  this  variation  is  significant, 
if  the  probe  orientation  is  known, 
the  variation  could  be  corrected 
mathematically. 

In  a  second  experiment, 
the  thermocouple  pro^  was 
placed  in  the  ultrasound  field  at 
normal  incidence  (with  the  probe 
axis  at  90°  to  the  wave 
propagation  vector.)  Each 
thermocouple  was  placed  in  the 


angle  between  axle  of  ineonation  and  probe  axle 
(degrees) 

Figure  34:  This  graph  shows  the  variation  in  rate  of  heating 
measured  by  three  different  thermocouples  within 
the  same  probe  as  the  angle  of  incidence  of  the 
ultrasound  field  was  varied  with  respect  to  the 
probe  axis. 
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focus  of  the  ultrasound  transducer  and  then  heated.  After  heating,  the  probe  was  rotated  about  its 
axis  and  the  experiment  was  repeated.  This  was  performed  for  each  of  the  seven  thermocouples 
in  the  probe.  The  results  of  this  experiment  are  shown  in  Figure  35.  Note  that  there  are 
significant  variations  with  angle  of  rotation  and  from  thermocouple  to  thermocouple.  This 
behavior  is  attributed  to  a  semi-random  placement  of  the  0.3mm  thermocouple  junction  within 
the  0.58mm  inside  diameter  of  the  tube.  In  order  to  better  understand  this  behavior,  the  finite 
difference  model  developed  previously  was  used  to  analyze  the  system. 

Finite  Difference  Model 


A  finite  difference  model  of  the  plastic  sheathed  thermocouple  probe  was  created  using  the 
forward  difference  technique  and  the  Matlab©  software  package.  The  model  geometry 
represents  a  thermocouple  probe  placed  on  the  axis  of  insonation  of  a  power  density  field 
generated  by  a  hyperthermia  treatment  transducer.  The  finite  difference  model  was  created  to 
simulate  the  transient  response  of  the  thermocouple  probe  at  the  onset  of  insonation.  Figure  36 
illustrates  the  2-D  cylindrical  coordinate  system  used  to  model  a  cross  section  of  the  probe  in  a 
water  bath.  Due  to  symmetiy  of  probe  and  incident  ultrasound  only  half  of  the  total  cross  section 
needed  to  be  modeled.  The  program  allowed  for  variation  of  the  model  radius  as  well  as  the 


number  of  nodes  in  the  radial 
direction.  The  number  of  angular 
divisions  could  also  be  varied  to 
efficiently  produce  meaningful 
results. 

Use  of  a  2-D  model  instead 
of  a  3-D  model  was  justified  because 
the  wide  energy  field  of  the  modeled 
treatment  transducer  produces  a 
radial  intensi^  gradient  near  the  axis 
of  insonation  that  is  not  steep.  The 
smooth  and  nearly  constant  intensity 
profile  proximal  to  the  cross  section 
of  interest,  which  is  centered  on  the 
axis  of  insonation,  and  the  low  heat 
conduction  coefficient  of  the  plastic 
produce  a  smooth  and  nearly 
constant  temperature  profile  along 
the  axis  of  the  probe  proximal  to  the 
cross  section  of  interest.  This 
smooth  and  nearly  constant 
temperature  profile  along  the  axis  of 
the  probe  near  the  cross  section  of 
interest  does  not  generate  extensive 
heat  conduction  down  the  probe 
relative  to  the  heat  produced  by  the 
ultrasound  at  the  cross  section  of 
interest.  The  nearly  constant 
intensity  within  a  few  millimeters  of 
the  axis  of  insonation  also  justifies  the  use  of  a  constant  incident  intensity  in  the  model. 

The  general  two-dimensional  cylindrical  heat-conduction  equation  is 


d^T  ^  1  d^T  1  d^T  q 
dr^  r  dr  ^  k 


a  dt 


Equation  7 


The  heat  generation  in  the  plastic  is  due  to  the  absorption  of  ultrasound  or 


^  ~  ^US^pt 


Equation  8 


where  aus  is  the  absorption  coefficient  of  the  plastic  sheathing  material  and  Ip,  is  the  intensity  of 
ultrasound  m  the  plastic.  However,  there  is  an  impedance  mismatch  between  the  plastic  and 
water  and  the  direction  of  ultrasound  propagation  is  not  perpendicular  to  the  surface  of  the 
plastic,  excluding  where  0=0.  Losses  of  ultrasound  intensity  at  the  water  plastic  interface  were 
accounted  for  using  the  angularly  dependent  reflection  coefficient  [Auld]  and  Equation  8 
becomes 
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Equation  9 


fV  .  -V 

/ cos(sin  (sin(^)  •  Cj  /c, ))  /cos(0) 

V  ,  +v 

/cos(sin"'(sin(^)c2/c,))  /cos{0) 


where  Z2  and  C2  are  the  impedance  and  speed  of  sound  in  plastic  respectively,  Z\  and  ci  are  the 
impedance  and  speed  of  sound  in  water  and  I  is  the  temporal  average  intensity  of  the  ultrasound 
incident  to  the  probe.  Intensity  loss  due  to  attenuation  in  the  plastic  as  the  ultrasound  propagated 
were  not  accounted  for  because  of  the  small  distance  actually  traveled  within  the  plastic  and  the 
inherent  estimation  of  intensity  in  the  plastic.  Estimation  of  the  impedance  of  the  plastic  and  of 
the  effect  of  total  reflection  of  ultrasoimd  intensity  at  the  plastic-air  interface  resulted  in 
imcertainty  of  the  actual  intensity  field  in  the  plastic. 


The  nodal  forward  difference  representation  of  the  general  2-D  cylindrical  heat-conduction 
equation  with  the  described  ultrasonic  absorption  is 


1  -  2a^ ,  A/| 


+  a^_,A/| 


TP  4.  TP 

Ar^ 


TP  _  TP 

^  ^  e,r-\ 

2rAr 


r^A0-  T/t, 


0,r 


where  each  node  carries  its  own  value  for  a,  rate  of  q,  and  k. 


Equation  10 


To  determine  an  appropriate  time  step  for  the  model  the  value  of  At  was  determined  after 
setting  the  following  quantity  to  zero: 


\-2agM 


Ar^ 


Equation 


This  value  of  At  was  the  largest  time  step  allowed  which  did  not  result  in  an  unstable  model.  To 
avoid  the  problem  of  dividing  by  zero  the  temperature  at  the  node  r  =  0  was  found  by  averaging 
all  of  the  surrounding  nodal  temperatures.  As  an  initial  condition  all  of  the  nodes  at  t  =  0  were 
set  to  room  temperature,  23  degrees  Celsius.  As  a  boundaiy  condition  the  temperature  at  the 
outer  radius  of  the  water,  R,  was  maintained  at  the  initial  condition  temperature.  Ultrasound 
intensity  was  modeled  as  a  step  function  beginning  at  t  =  0. 

Rcsiilts 


Figure  37  shows  the  heating  characteristics  of  the  finite  difference  model  of  the  insonated  clinical 
thermocouple  probe  after  one  second  of  ultrasound  insonation.  It  can  be  seen  that  the  quadrant  of 
the  probe  which  ‘  sees”  the  ultrasound  intensity  is  the  portion  of  the  probe  that  heats  due  to 
ultrasound  absorption.  The  heating  of  the  far  side  of  the  sheath  and  heating  of  the  water  are  due 
to  conduction.  It  can  be  seen  from  the  figure  that  the  temperature  profile  from  the  insonated  side 
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Figure  38:  This  shows  the  predicted  variation  in  temperature  for  points  at 

various  distances  from  the  central  axis  of  the  probe  as  the  probe  is 
rotated  about  its  axis. 


•  A  series  of  new  probes  were  designed,  constructed  and  tested. 

•  A  flnal  design  was  proposed  for  future  work. 

A  series  of  probes  have  been  built  focussing  on  forcing  the  thermocouples  into  the  center 
of  the  probe.  These  have  reduced  the  variation  somewhat  and  made  the  variation  more  reliably 
sinusoidal,  but  a  transducer  has  not  been  built  to  date  that  reduces  this  variation  below  +/- 10%. 
A  new  design  has  been  proposed  that  is  intended  to  further  reduce  the  variation.  [29] 


Subsystem  3;  Model  Development 

The  model  was  completed  in  the  3"*  year.  No  work  was  accomplished  on  this 
subsystem  during  the  4**"  year. 


due  to  ultrasound  absorption.  The  heating  of  the  far  side  of  the  sheath  and  heating  of  the  water 
are  due  to  conduction.  It  can  be  seen  from  the  figure  that  the  temperature  profile  from  the 
insonated  side  of  the  sheath  to  the  far  side  is  decreasing  significantly.  This  supports  the 
postulation  that  the  transient  response  of  a  sensor  within  the  clinical  thermocouple  probe  is 
largely  dependent  on  the  position  of  the  sensor  within  the  probe. 

According  to  the  finite  difference  heating  pattern  of  the  insonated  clinical  thermocouple 
probe  it  would  be  expected  that  if  the  probe  were  rotated  about  its  own  axis  that  the  effective 
absorption  parameter  would  be  highly  dependent  on  rotation.  It  would  be  expected  that  a  plot  of 
the  effective  absorption  parameter  magnitude  versus  the  respective  rotation  for  an  off  centered 
sensor  would  be  sinusoid-like  in  nature. 

Figure  38  is  a  plot  of  the  normalized  temperature  at  distances  from  the  probe  axis  versus 
theta  after  one  second  of  insonation.  It  is  clear  that  the  final  temperature  at  a  distance  from  the 
axis  has  a  sinusoid-like  dependence  with  theta,  and  increasing  the  off  center  distance  also 
increases  the  magnitude  of  oscillation.  Because  the  effective  absorption  parameter  is  determined 
by  a  linear  curve  fit  of  temperature  rise  the  variations  in  temperature  rise  are  directly 
proportional  to  the  determined  effective  absorption  parameter.  Consequently,  Figure  38  shows 
that  a  probe  with  a  sensor  off  centered  by  0. 1-mm  will  result  in  ±  22%  variation  of  the  effective 
absorption  parameter  with  probe  rotation  and  a  sensor  off  centered  by  0.2-nim  will  result  in 
approximately  ±43%  variation  of  the  effective  absorption  parameter  with  probe  rotation. 


Subsystem  4;  Inverse  Tcchniouc 

Subsystem  integration  completed. 

•  Technique  used  successfully  to  optimize  the  input  power  model. 

During  the  fourth  year,  the  inverse  optimization  technique  discussed  in  previous  reports 
was  integrated  with  the  other  subsystems,  “^is  provides  the  anticipated  capability  to  adjust  the 
tissue  attenuation  values  to  minimize  the  error  between  predicted  power  deposition  values  and 
experimentally  measured  values.  This  capability  was  not  implemented  however.  The  ultrasonic 
attenuation  measurement  was  shown  to  measure  attenuation  with  an  error  of +/-  10%.  These 
measurements  are  the  input  nodal  attenuation  values  for  the  model  (before  optimization).  Since 
the  experimental  SAR  measurements  were  not  capable  of  an  error  as  small  as  this,  the 
optimization  would  provide  no  benefit.  The  capability  for  optimization  is  available  when  the 
SAR  measurement  technique  improves  sufficiently  to  provide  a  benefit. 

The  optimization  routines  were  used  for  an  alternative  application  that  does  provide  an 
increase  in  accuracy  of  the  model.  This  involves  the  optimization  of  the  model  prediction  for  the 
ultrasound  field  produced  by  the  treatment  transducer.  This  application  is  described  briefly 
below  (and  more  completely  in  [30]): 

Optimizing  the  Power  Input  Field 

Numerical  predictions  of  the  ultrasound  field  produced  by  a  focused  transducer  show 
small  but  significant  variations  from  the  results  of  experimental  measurements  when  the 
transducer  is  modeled  as  a  piston.  When  the  normal  velocities  on  the  surface  of  the  modeled 
transducer  are  allowed  to  vaiy,  a  more  accurate  model  results.  In  this  work  the  acoustic  field  of  a 
500  kHz  spherically  focussed  transducer  was  experimentally  mapped  using  a  3  axis  motion 
device  and  a  hydrophone.  A  model  of  the  transducer  was  developed  using  a  uniform  distribution 
of  point  sources  on  the  surface  of  spherical  section.  The  particle  velocities  at  each  point  were 
defined  based  on  control  points  for  a  1  -D  cubic  spline  (radial  symmetry  is  assumed).  By  varying 
the  number  and  position  of  control  points,  different  velocity  functions  for  the  transducer  face  can 
be  created  (A  single  control  point  produces  piston-like  behavior,  two  points  result  in  a  linear 
variation,  etc.).  The  control  point  velocities  were  optimized  using  a  conjugate  gradient 
minimization  algorithm  (Polak-Ribiere).  A  piston  model  produces  results  similar  to  the  existing 
USSPAR  model.  Using  four  or  more  control  points  reduces  the  error  between  model  and  data  by 
ten  percent.  The  resulting  velocity  function  is  reminiscent  of  a  Bessel  function,  perhaps 
suggesting  that  the  transducer  has  a  secondaiy  drum-like  vibrational  mode  in  addition  to  the 
expected  piezo-electric  mode.  Further  improvements  may  be  obtained  by  removing  the 
symmetry  constraint  using  a  two-dimensional  spline  fit. 

Introduction 

Focused  ultrasound  transducers  are  used  in  hyperthermia  cancer  treatment  therapy.  This 
treatment  when  combined  with  radiation  has  been  shown  to  nearly  double  the  response  rates  over 
radiation  alone  [31].  The  goal  of  the  treatment  is  to  uniformly  heat  a  tumor  to  43°  C  [32]  while 
minimizing  the  temperature  rise  in  the  surrounding  tissue  (Figure  23).  To  address  this  problem, 
numerical  models  have  been  created  to  simulate  the  hyperAermia  temperature  field  [33-37]. 

Since  such  models  require  an  estimate  of  the  power  absorbed  by  the  tissue,  work  is  being  done  to 
create  a  more  accurate  power  model.  The  calibration  of  the  treatment  transducer  is  the  first  step 
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in  reliably  predicting  the  absorbed  power  field.  One  method  for  doing  this  is  explained  by 
Swindell  [38,39].  His  model  (USSPAR),  although  it  is  fast,  is  based  on  the  assumption  that  the 
velocities  on  the  transducer  face  are  constant.  More  recently,  however,  models  have  been 
developed  which  allow  for  variable  surface  velocities.  One  uses  superimposed  Gaussian  beams 
to  represent  the  velocity  profile  [40],  and  others  break  the  transducer  into  elemental  areas 
[41,42]. 

In  this  work,  the  transducer  is  represented  as  a  collection  of  point  sources.  The  velocities  at 
each  point  are  defined  based  on  control  points  for  a  1-D  cubic  spline  assuming  radial  symmetry. 
By  varying  the  number  of  control  points,  different  velocity  functions  can  be  created  (A  single 
control  point  produces  piston-like  behavior,  two  points  result  in  a  linear  variation,  etc.).  This 
method  was  employed  to  simplify  the  inverse  problem  so  that  large  focussed  treatment 
tran^ucers  could  be  calibrated  more  quickly  (Figure  39).  A  conjugate  gradient  minimization 
routine  was  used  to  fit  the  model  to  data  generated  by  USSPAR  as  well  as  measured  data. 
Experimental  Setup 

Data  was  measured  using  a  hydrophone  and  the  setup  shown  in  Figure  40.  The  500  kHz. 
transducer  shown  in  Figure  39  was  mounted  to  a  3  axis  motion  control  device  which  was 
controlled  using  a  PC.  The  hydrophone  was  mounted  to  a  fixed  fiame  inside  a  tank  filled  with 
degassed  water.  Output  fi-om  the  hydrophone  was  amplified  and  read  from  an  oscilloscope.  The 
voltage  amplitude  was  measured  and  recorded  at  each  data  point  along  with  the  corresponding 
scale  factor.  The  intensity  was  calculated  using  the  experimentally  determined  calibration  factor 
for  the  hydrophone  at  this  frequency  (500  kHz.):  4  =  1 .595  •  where  is  the  peak  to  peak 
voltage  measured. 


Cross  Section  of  Transducer 


7cm 


Figure  39.  This  shows  the  construction  of  a  typical  hyperthermia  treatment 
transducer.  The  crystal  is  suspended  by  it’s  edges  and  air  backed 
to  maximize  the  energy  transmitted  in  the  preferred  direction. 
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Theory 

The  Transducer  Model 

The  transducer  was  modeled  as  a  group  of  point  sources  uniformly  distributed  on  a  section  of  a 
sphere.  This  was  accomplished  by  mapping  a  uniform  set  of  points  in  a  plane  onto  a  sphere. 
Since  the  transducer’s  radius  of  curvature  is  much  larger  than  it’s  radius,  the  error  in  this 
approximation  is  small.  The  distance  between  point  sources  was  limited  to  be  about  a  half  a 
wavelength.  The  pressure  field  caused  by  the  transducer  was  calculated  by  differentiating  the 
Rayleigh-Sommerfeld  diffraction  integral  producing  the  following  equation; 


where  ^  =  <y/c  =  2;r//l  is  the  wave  number,  2  is  the  wavelength,  Z  =  pc  the  acoustic 
impedance  of  the  propagating  medium,  v(r^.)=  VjC*"*  is  the  normal  particle  velocity  on  ds,  and 

^ij  “  h  “  oi  distance  from  a  point  in  the  field,  r  =  r^  to  a  point  on  the  transducer,  r  =  rj.  If 
ds  =  {5Xf  is  constant  for  N  point  sources  on  S  where  SX  «  r  (not  close  to  transducer)  the 

pressure  can  be  approximated  by:  =  iZS  X^Vj - where  ^  is  the  number  of 

7=1 

wavelengths  between  points.  If  ^  ^r^jX,  this  simplifies  to: 
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N  'V  P.P. 

Pi  =  iZS^^Vj—z —  =  iZS^^Ef/Vj  where  Ey  =  — .  Since/^  =  -^(no  sum  on  i),  we  get 

M  J=i  %  ^ 

N  N  _  _  _ 

i,=zs‘y;zv,e,  EyVj  (Pi  and  are  the  complex  conjugates  to  and  Ey.  respectively).  To 

^=I  >=1 


reduce  computation  time  (by  decreasing  the  dimension  of  parameter  space),  and  to 
simultaneously  enforce  a  relatively  smooth  velocity  distribution,  a  1-D  cubic  spline  routine  is 
used  to  interpolate  velocities  at  each  point  on  the  transducer  (radial  symmetry  assumption).  The 
velocity  profile  on  the  transducer  face  can  then  be  modified  by  changing  the  velocities  at  the 
control  points.  The  cubic  spline  routine  generates  a  transformation  matrix  which  maps  model 
parameters  (spline  control  point  velocities)  to  transducer  point  velocities.  If  this  matrix  ( Tj,. )  is 


M 

defined  by:  v^.  =  ,  the  intensity  is  then: 

t=i 

N 

tensors  are  defined  such  that  ,  and 

>=i 


MM  (  N  N  _  ^ 

1, = zs-'Z'Z'x,  If 


A,n  s  ZS^ByBy^  (no  sum  on  /),  the  intensity  is 


M  M 

then  just  7,  'Si;  m,A,y^m,^ .  When  iterating,  only  the  model  parameters  (m )  will  change  so 


/=1  i=l 


[M^DxM] 

the  third  order  tensor  (  A,y^  )  needs  to  be  calculated  only  once  and  can  be  easily  stored  in 
memory  because  M  is  small. 


Model  Optimization  Technique 


The  Conjugate  gradient  method  is  used  to  minimize  the  error  function,  <1>  s  -^l/  - 
where  is  the  measured  intensity  data.  Such  methods  have  been  shown  to  produce  favorable 
results  for  a  wide  range  of  applications  [13-16].  The  algorithm  used  here  is  based  on  a 
procedure  in  Numerical  Recipes  w  C  [17].  The  search  direction  is  defined  to  be 

Cj?'  —  • 

■^/+i  =  “g/+i  +  ySi  where  y  =  — (Polak-Ribiere).  The  search  starts  in  the  steepest 

g,  Si 


descent  direction  (go  =  g, )  and  then  each  subsequent  direction  is  conjugate  to  all  previous  ones. 
The  gradient  is  found  by  differentiating  with  respect  to  the  model  parameters,  m: 

T  do  r  dl 

S  =  ^  =  (^  -  4)  •  The  Jacobian  is  simply  J  =  —  =  2mA .  The  step  length,  a ,  for  a  1-D 
am  dm 

line  minimization  is  found  by  using  a  linear  approximation  for  the  Intensity 
( I (m  +  as)»  I(m)+ aJs ).  The  error  function  becomes  <I>  =  -^l/ — 7^  +  aJs^ .  Minimizing  with 

do  T 

respect  to  step  length  gives  —  =  (7  -  7^  +  aJsf  Js  =  0 .  The  step  length  is  then  a  =  . 
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Results 


To  initially  verify  the  accuracy  of  the  model,  theoretical  data  was  generated  using 
USSPAR,  an  existing  power  model  which  calculates  the  field  under  the  assumption  that  the 
velocities  on  the  surface  of  the  transducer  are  constant  [38,39], 

The  piston  velocity  profile  of  both  USSPAR  and  the  piston  model  with  one  spline  point 
were  almost  identical.  Thus,  the  piston  model  fit  USSPAR  well 

(misfit,  =  ||/  -4||/||4||  *  0.0022 )  but  neither  did  a  good  job  of  fitting  the  data  (A^  *  0. 1743 ). 
When  four  spline  points  were  used,  however,  the  fit  improved  considerably  {Mf  »  0.0634 ). 
Beyond  four  points,  the  misfit  didn’t  decrease  much  (The  12-point  model  had  Mf  «  0.0620 .), 
but  the  curve  became  smoother  and  consequently  more  believable.  As  the  number  of  points  and 
point  locations  changed,  the  general  curve  shape  remained  the  same.  Finally,  reducing  the 
transducer  point  spacing  (<y  «  0.4 )  failed  to  improve  the  fit  ( A^  »  0.0622 ),  indicating  that  the 

model  points  were  sufficiently  dense  to  produce  an  accurate  representation  of  a  continuous 
surface. 


Trarcduoer  VflIocHy  Profile  [Pfeion  Model) 


Ffl  lor  Pisfon  Model  fTo  USSPAR  &  Dala) 


DiBfarvn  From  Gerrier  of  Field,  X  (cm) 


Figure  4 1 .  When  the  model  uses  a  piston-like  velocity  profile  it  behaves  identically  to  the 
USSPAR  model,  providing  a  relatively  poor  match  to  the  experimental  data. 
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Figure  42:  When  the  surface  velocity  profile  on  the  transducer  face  is  allowed  to  vary,  the 
misfit  between  model  data  and  experimental  data  is  reduced. 


Discussion 

The  main  goal  of  this  research  was  to  develop  an  accurate  model  of  a  focussed 
transducer,  but  it  was  also  hoped  to  obtain  a  realistic  estimate  of  the  transducer’s  actual  surface 
velocity  field.  In  searching  for  the  actual  field,  one  has  the  advantage  of  using  physiologic 
feasibility  to  help  define  the  solution  [48],  Also,  when  the  model  uses  a  realistic  velocity  profile, 
it  is  more  likely  to  predict  realistic  values  outside  the  range  of  the  data  used  in  the  optimization 
procedure.  One  of  the  difficulties  in  optimization  is  knowing  when  to  quit.  If  the  data  is  fit  too 
closely,  the  model  may  be  fitting  experimental  noise  and  consequently  give  poor  results.  For 
these  reasons,  a  smooth  solution  was  one  of  the  goals  (It  is  expected  that  large  variations  in 
velocities  would  tear  the  transducer  apart).  It  is  also  a  good  sign  when  varied  initial  conditions 
converge  to  a  similar  solution.  For  instance,  when  more  spline  points  are  used,  the  solution 
should  not  change  drastically  but  should  be  refined.  Originally,  the  resulting  curves  varied  quite 
dramatically  when  more  points  were  included.  It  was  then  discovered  that  the  convergence 
tolerance  was  too  small.  Once  the  routine  was  modified  to  quit  when  the  misfit  began  to  level 
out,  the  solutions  all  began  to  look  alike  regardless  of  the  number  of  points  used. 


^8 


TiarBduoer  VGlocity  Pfolile  (12-^1.  Model) 
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Figure  43;  Using  12  spline  points  provides  a  significantly  improved  fit. 

The  refined  velocity  function  appears  to  be  reminiscent  of  a  Bessel  function,  which  may 
suggest  that  the  transducer  has  a  secondary  drum-like  vibrational  mode  in  addition  to  the 
expected  piezo-electric  mode.  This  may  be  due  somewhat  to  the  construction  of  the  air-backed 
transducer  used  (Figure  39).  The  relatively  unconstrained  moimting  setup  should  permit  some 
resonant  behavior. 

When  reducing  the  transducer  point  spacing  failed  to  improve  results,  it  was  decided  that 
the  remaining  error  must  lie  mainly  in  the  measured  data  rather  than  in  the  model.  The  model’s 
main  assumption  is  that  the  transducer  can  be  modeled  as  a  collection  of  points  if  they  are 
uniformly  spaced  and  close  together.  Consequently,  increasing  the  number  of  points  should 
reduce  the  error  in  the  model. 

Conclusions  and  Recommendations 

Although  it  is  still  uncertain  whether  the  model’s  surface  velocity  field  agrees  with 
reality,  there  is  reason  to  believe  that  it  is  a  good  approximation  to  reality  since  it  is  smooth, 
consistent  and  fits  the  data  quite  well.  It  has  also  been  shown  that  allowing  the  surface  velocities 
to  vary  can  be  beneficial.  This  particular  model  is  not  expected  to  be  accurate  at  close  range, 
however.  It  is  suited  particularly  for  large  focussed  treatment  ultrasound  transducers  where 
accuracy  in  the  veiy  near  field  is  not  required.  This  method  is  quite  simple  to  apply  and  not 
computationally  expensive.  All  the  code  was  written  using  Matlab  5.0  and  the  12-point  model 
(N=1725,  D=55)  took  less  than  30  seconds  to  converge  on  a  Sun  SPARC  20  workstation. 


r 


Further  improvements  may  be  obtained  by  removing  the  symmetry  constraint  using  a  two- 
dimensional  spline  fit.  A  more  detailed  discussion  of  the  optimization  results  can  be  foimd  in  []. 


•  Given  the  level  of  scatter  in  the  experimental  SAR  measurements,  the 
optimization  model  was  not  used  to  adjust  nodal  attenuation  values. 

Conclusions  and  Recommendations 

The  results  of  this  project  represent  a  significant  step  forward  in  accuracy  for  a  practical, 
clinical  ultrasound  power  deposition  model  for  hyperthermia  treatment  plaiming.  While  this  is 
true,  accurate  power  deposition  is  only  one  piece  of  a  complete  model.  This  model  must  be 
integrated  into  a  complete  thermal  model.  Other  issues  will  continue  to  limit  the  accuracy  of  the 
thermal  models,  including  unknown  perfusion  values  and  changes  in  the  model  parameters  that 
result  from  the  treatment.  The  statement  of  work  included  anticipated  clinical  trials  of  the 
system  which  were  not  carried  out.  This  is  a  disappointment  and  results  from  the  dependence  of 
this  work  on  the  clinical  data  of  other  researchers.  Anticipated  clinical  tests  in  the  Hyperthermia 
Lab  at  the  University  of  Utah  Hospital  have  not  yet  been  conducted. 

The  system  developed  here  is  currently  being  implemented  on  the  hyperthermia 
treatment  facility  at  the  University  of  Utah  Hospital  and  clinical  trials  are  anticipated  in  future 
work.  The  implementation  of  the  system  developed  here  in  the  Hyperthermia  Lab  involves  two 
steps:  First,  the  hardware  must  be  installed  on  the  clinical  system.  This  was  completed  as  part 
of  the  final  testing  of  the  system.  Second,  the  model  integration  mentioned  above  must  be 
completed.  This  is  currently  being  investigated  by  researchers  in  the  Hyperthermia  Lab. 

As  a  final  recommendation,  it  is  anticipated  that  the  accuracy  and  value  of  the  model 
developed  here  can  be  improved  by  further  development  of  an  accurate  experimental  method  to 
measure  SAR. 
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